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Abstract 

This paper deals with the dead-water phenomenon, which occurs when a ship sails in a 
stratified fluid, and experiences an important drag due to waves below the surface. More gen- 
erally, we study the generation of internal waves by a disturbance moving at constant speed 
on top of two layers of fluids of different densities. Starting from the full Euler equations, 
we present several nonlinear asymptotic models, in the long wave regime. These models are 
rigorously justified by consistency or convergence results. A careful theoretical and numerical 
analysis is then provided, in order to predict the behavior of the flow and in which situations 
the dead-water effect appears. 

1 Introduction 

The so-called "dead-water" phenomenon has been first reported by Fridtjof Nansen [3^, describing 
a severe (and then inexplicable) decrease of velocity encountered by a ship, sailing on calm seas. 
Bjerknes, and then Ekman |25j soon explained that this phenomenon occurs as the ship, on 
top of density-stratified fluids (due to variations in temperature or salinity), concedes a large 
amount of energy by generating internal waves. Our work is motivated by the recent paper of 
Vasseur, Mercier and Dauxois |42j . who performed experiments that revealed new insights on the 
phenomenon. 

Relatively little consideration has been given to this problem, after the early works of Lamb [32] 
and Sretenskii 54 . Miloh, Tulin and Zilman in |44| produce a model and numerical simulations 
for the case of a semi-submersible moving steadily on the free surface. The authors assume that 
the density difference between the two layers is small, an assumption that is removed by Nguyen 
and Yeung [33]. Motygin and Kuznetsov [3S] offer a rigorous mathematical treatment of the 
issue when the body is totally submerged in one of the two layers, and Ten and Kashiwagi [55] 
present a comparison between theory and experiments, in the case of a body intersecting the 
interface as well as the surface. Finally, we would like to cite Lu and Chen [33] for a more general 
treatment of the problem. All of these works use linearized boundary conditions, that rely on the 
assumption that the amplitudes of the generated waves are small. The linear theory is convenient 
as it allows to obtain the flow field by a simple superposition of Green functions, replacing the 
moving body by a sum of singularities. The integral representation of the flow, as well as the 
wave resistance experienced by the body, are therefore found explicitly. However, the smallness 
assumption on the wave amplitudes, that is necessary to the linear theory, is quite restrictive, and 
the experiences in [42] clearly exhibit nonlinear features. Our aim is to produce simple nonlinear 
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models, that are able to predict the apparition and the magnitude of the dead-water effect, and 
recover duly most of its key aspects. 

In this paper, we introduce several different asymptotic models, corresponding to two different 
regimes (size of the parameters of the system). Each time, we assume that the depth of the fluid 
layers are small when compared with the wave lengths (shallow water). In the one-layer case, the 
equations obtained with this single assumption are the shallow-water or Saint Venant equations 
(at first order), or the Green-Naghdi extension (one order further, therefore involving nonlinear 
dispersion terms). See Wu, Chen [57j and references therein for a numerical treatment of the waves 
generated by a moving ship in the one layer case, using the Green-Naghdi equations. The two 
regimes we study carry additional smallness assumptions, that allow to substantially simplify the 
models. The first regime considers the case of small surface deformations, and small differences 
between the densities of the two layers. Such additional assumptions are very natural in the 
oceanographic framework, and have been frequently used in the literature (see [121 [Ml US] for 
example). In the second regime, we assume that the magnitude of the produced internal waves, 
when compared with the depth of the two layers, are small and of the same order of magnitude 
as the shallowness parameter. This regime, known as the Boussinesq regime, is particularly 
interesting as it allows models with competing dispersion and nonlinearity. Along with the coupled 
Boussinesq-type model, we introduce the KdV approximation, which consists in a decomposition 
of the flow into two waves, each one being a solution of an independent forced Korteweg-de 
Vries equation (fKdV) . The fKdV equation has been extensively studied in the framework of the 
one-layer water wave problem (where a moving topography, or pressure, is the forcing term that 
generates waves); see [SSI UU [37] , for example. 

The system we study consists in two layers of immiscible, homogeneous, ideal, incompressible 
fluids only under the influence of gravity. The bottom is assumed to be flat, and we use the 
rigid-lid approximation at the surface. However, the surface is not flat, but a given function, 
moving at constant speed, that reflects the presence of the ship. Moreover, as we are interested 
in unidirectional equations (the fKdV equations), we focus on the two-dimensional case, i.e. 
the horizontal dimension d = 1. However, our method could easily be extended to the three- 
dimensional configuration. Starting from the governing equations of our problem, the so-called 
full Euler system, and armed with the analysis of [231 [23 in the free surface case, we are able to 
deduce asymptotic models for each of the regimes presented above. Each of the models presented 
here are justified by a consistency result, or a convergence theorem (in a sense precisely given in 
Section 12.31 page II Op . We compute numerically the fully nonlinear system of the first regime, as 
well as the KdV approximation, which allows us to investigate the effect of different parameters 
of the system, such as the velocity of the boat, or the depth ratio of the two layers. The wave 
resistance encountered by the ship is also computed, so that we are able to predict in which 
situations the dead-water effect shows up. 

Our models allow to recover most of the known aspects of the dead-water phenomenon (see [HI 
mi [in]), especially: 

i. transverse internal waves are generated at the rear of a body while moving at the surface; 

ii. the body suffers from a positive drag when an internal elevation wave is located at its stern; 

iii. this effect can be strong (as the generated wave reach large amplitudes) near critical Froude 
numbers, that is when the velocity of the body approaches maximum internal wave speed. 
The effect is always small away from the critical value; 

iv. the maximum peak of the drag is reached at slightly subcritical values. 
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We would like to emphasize that these considerations are not new. In particular, the fact that the 
drag experienced by the ship is strong only when its velocity approaches the maximum internal 
wave speed has been observed through experiments [251 142] , and the maximum peak of the drag 
being reached at slightly subcritical values is predicted by linear models [HI US]- The main 
contribution of our work is to offer simple adapted models, that are rigorously justified for a 
wide range of parameters (in particular, in contrast with linear models, our models allow large 
amplitude waves, in agreement with experiments). 

Our models allow to easily investigate the dead-water phenomenon for several values of the 
parameter. For example, we numerically compute the system for several values of the depth 
ratio between the two layers. It appears that, in disagreement with the intuition, the dead-water 
phenomenon is stronger when the upper layer is thicker than the lower one. This observation is 
new, as far as we know. This curiosity is due to the fact that the magnitude of the drag suffered by 
the body does not depend on the distance between the body and the generated internal wave, but 
rather on the amplitude and sharpness of the latter. A thicker upper layer allows the generation 
of internal elevation waves reaching very large amplitudes, when compared with the size of the 
body, which strengthen the dead-water phenomenon. 

The behavior of the system, depending on the depth ratio of the two layers of fluid and the 
normalized velocity of the body, is summarized in Table [TJ below. 

However, one peculiar phenomenon, that is described in details in |42j . is not recovered by our 
models. Indeed, the dead-water effects exhibits a somewhat periodic behavior, where during each 
period, a wave is generated, slows down the ship, and then breaks. This discrepancy between our 
simulations and the experiments is due to the fact that our models are based on the assumption of 
a constant velocity for the traveling body, while their experiments are conducted with a constant 
force brought to the body (5S1I12]. The latter setting is of course more natural, but the constant 
velocity hypothesis is crucial in our analysis (and is, in fact, consistently assumed within all 
existing models in the literature |25l 1441 [55] ). We perform a numeric experiment which roughly 
matches this setting (the velocity is then adjusted at each time step as a function of the drag 
suffered by the ship), and recover the oscillating behavior. This point is discussed in more details 
on page [THl 

Outline of the paper : In Section [531 we introduce the governing equations of our problem: 
the so-called full Euler system. A specific care is given to an adapted rescaling, that leads to 
the dimensionless version of the full Euler system (|2.4p . The models presented here are asymp- 
totic approximations of this system. In Section 12.31 we present precisely the different means 
of justification of our models, as well as the regimes at stake. The main tool that we use in 
order to construct our models, as well as a broad, strongly nonlinear model ([2. 71) are presented 
in Section [231 

The simpler models we study are deduced from ([2.7[) . using the additional assumptions of 
the regimes considered. Strongly nonlinear models (13. ip and (13. 3p are introduced in Section [31 
and justified with consistency results. Numerical simulations are then displayed and discussed 
in Section [3.21 The weakly nonlinear models, i.e. the Boussinesq-type system (|4.2p (and its 
symmetrized version (14. 3p ). and the KdV approximation (j4.12p . are presented in Sections 14.11 
and 14.21 respectively. The convergence of their solutions towards the solutions of the full Euler 
system are given in Propositions 14.11 and 14.21 respectively. An in-depth analysis of the forced 
Korteweg-de Vries equation, and its consequences to the dead-water effect, is then displayed in 
Section [131 

Finally, in Section [SI we conclude with an overview of the results presented here, together 
with suggestions for future work. 

Calculations that lead to (|2.7p are postponed until Appendix [XI and Appendix [Bl contains the 
proof of Proposition 14.21 Finally, Appendix [O is devoted to the analysis of the wave resistance 



Velocity of 
the ship 


subcritical case 


critical case 


supercritical case 


Depth ratio 
between 
the layers 




thicker lower layer 


thicker upper layer 




Regime 1 


The generated waves are 
small, and conclusions of 
Regime 2 apply. 


Generation of dull 
elevation-depression wave 
below the body. 

Moderate wave resistance. 

See Figure [21 pagelTTI 


Generation of a sharp 
elevation wave below the 
body. 

Strong wave resistance. 
See Figure m pagelTTI 


The generated waves are 
small, and conclusions of 
Regime 2 apply. 


Regime 2 


Continuous generation of 
small up-stream 
propagating and 
down-stream propagating 
waves. 

Very weak wave 
resistance. 


Periodic generation of 
up-stream propagating 
depression waves with 
very large time-period0 

Oscillating wave 
resistance, with positive 
mean. 


Periodic generation of 
up-stream propagating 
elevation waves with very 
large time-period." 

Oscillating wave 
resistance, with mean 
zero. 


Generation of small 
down-stream propagating 
waves. Convergence 
towards a steady state. 

No lasting wave 
resistance. 




See Figure [6l page [28l 


See Figure [9l pagel30l 


See Figure m page [30] 


See Figure III page|28| 



Table 1: Behavior of the flow in the two studied regimes, depending on the velocity of the body and the depth ratio of the fluids 

"As discussed in Section 14.3.21 the time period of the generation of the waves, and therefore of the oscillation of the wave resistance, may be larger than the 
time-range of validity of the model. 
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encountered by the ship, and the numerical schemes used in the simulations are presented and 
justified in Appendix [Dl 

Notations : We denote by C(Ai, A2, . . .) any positive constant, depending on the parameters 
Ai, A2, . . ., and whose dependence on Xj is always assumed to be nondecreasing. 

Let f{xi, . . . , Xd) be a function defined on M'^. We denote by dx^f the derivative with respect 
to the variable Xi. If / depends only on Xi, then we use the notation 

dxi 

We denote by = the Hilbert space of Lebesgue measurable functions F : K — >■ 

1 /2 

with the standard norm \F\^ — ( Jj^|_F(a:)p da; ) < ck). Its inner product is denoted by 

Then, we denote by H'^ ~ H'^{M.) the L^-based Sobolev spaces. Their norm is written 
I • = [a* • with the pseudo-differential operator A = (1 — 9^)^/^. It is also convenient 
to introduce the following norm on the Sobolev space H^~^^ (for e > a small parameter): 

Let < r < 00 and f{t, x) be a function defined on [0, T] xK. Then one has / G L°°{[0, T]-H') 
if / is bounded in H", uniformly with respect to f € [0,r]. Moreover, / € W''-'°°{[0,T]; H") if 
/, dtf e L°°{[0,T];H'^). Their respective norm is written | • and | • ■ 



2 Construction of asymptotic models 

Our framework, the formulation of the system using Dirichlet-Neumann operators, as well as the 
techniques used for constructing the asymptotic models, are not new. We follow here the strategy 
initiated in [51 [SJ [7] in the one-layer (water wave) configuration, and extended to the bi-fluidic 
case in [HlllSlUS]- The major difference here is that the surface is not free, but rather a non-flat 
rigid lid, moving at constant velocity. 

We present in the following, as briefly as possible, the governing equations of our problem. The 
full details of the construction can be found in the references above. Our models are then derived 
starting from the non-dimensionalized version of these equations (the full Euler system (|2.4p ; see 
Section l2.2p . The idea is to use smallness assumptions on some parameters, that characterize 
natural quantities of the system, to construct the asymptotic models. Together with the means 
to justify our models, the regimes under study are precisely described in Section 12.31 Finally, 
the main tool that we use in order to construct our models, that is the asymptotic expansion 
of the Dirichlet-Neumann operators in terms of small shallowness parameter /i, is presented in 
Section A broad, nonlinear model (|2.7p is then deduced, and the precise systems we study 
are derived from it, using the additional assumptions of Regime [T] and Regime [U respectively in 
Section [3] and Section [H 



2.1 The basic equations 

The system we study consists in two layers of immiscible, homogeneous, ideal, incompressible 
fluids only under the influence of gravity (see Figure [T|) . Since we are interested in unidirectional 
models (the KdV approximation), we focus on the two-dimensional case, i.e. the horizontal 
dimension d = I. However, most of this study could easily be extended to the case d = 2, using 
the techniques presented here. 
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X 

Figure 1: Sketch of the domain 



We denote by pi and p2 the density of, respectively, the upper and the lower fluid. Since we 
assume that each fluid is incompressible, pi and p2 are constant and the velocity potentials (pi 
{i = 1,2), respectively associated to the upper and lower fluid layers, satisfy the Laplace equation 

dl(i>i + d% = 0. 

Moreover, it is presumed that the surface and the interface are given as the graph of functions 
(respectively, (^i{t,x) and C2(i, s;)) which express the deviation from their rest position (respec- 
tively, (a;, di) and {x, 0)) at the spatial coordinate x and at time t. The bottom is assumed to be 
flat, and the surface is flat away from the location of the ship moving at constant speed Cg, so 
that Ci matches the submerged part of the ship and is given by 

Cl{t,x) = (l{x-Cst). 

Therefore, at each time t > 0, the domains of the upper and lower fluid (denoted, respectively, 
0,\ and are given by 

n\ = { {x,z) eR'^ xR, C2{t,x) < z < di+Ci{x-Cst)}, 
nl = { (x, z) e R'' X R, -d2 < z < C2{t,x)}. 

The fluids being ideal, they satisfy the Eulcr equation (or Bernoulli equation when written in 
terms of the velocity potentials), and kinematic boundary conditions are given through the as- 
sumption that no fluid particle crosses the surface, the bottom or the interface. Finally, the set of 
equations is closed by the continuity of the stress tensor at the interface, which takes into account 
the effects of surface tension. 
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(2.1) 



Altogether, the governing equations of our problem are given by the following 



920, + 920, = 



in nl, i = 1,2, 

on Ti = {{x,z),z = di + Ciit,x)}, 



dtC.2 Vl + l^2^C2pa„2 0i = + |52;C2p9„202 on r2 = {(a;,z),z = (2(^,2;)}, 

on Tfe = {(x, z), z = -^2}, 

9xC2 \ 



9,02 = 



,v/i + |9.C2lV 

where 9„. is the upward normal derivative at 



on r2 



with 



5 we denote a > the surface tension coefficient and |P| the jump of the pressure at the interface: 

[P(i,a;)l = Ihn (P(t,a;, (2(^,2;) +e) - P(t, a;, C2(i, x) - e) ) . 

This system can be reduced into evolution equations located at the surface and at the interface. 
Indeed, when we define the trace of the potentials at the surface and at the interface 



V'i(t, = 0i(t, +Ci(i,x)), 



and 



-02 (i, a;) = 02(i,x,C2(t,a;)), 



10 then 01 and 02 are uniquely given by the Laplace equation, and the (Dirichlet or Neumann) 
boundary conditions (see Figure [2]) : 



( 92 + 92 ) 02 =0 in r!2, 
(2.2) { 02=^2 onF2, 

9,02 = on Ffc, 



( 92 + 92 ) 01 = in rii, 
and i 01 = ^1 on Fi, 

9„2 0i = 9„2 02 on F2. 




Figure 2: Laplace problems in the two domains. 



Therefore, the following operators are well-defined: 

Gi[Ci,C2] (^1,^-2) = x/l + |9,Ci|2(9„,0i) 

G2[C2]V'2 = Vl + |5.C2|2(9„,02)|„,^ , 

^^[Ci,C2](V'i,V'2) ^ ^.(0iL.c. )• 
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Then, the chain rule and straightforward combinations of the equations of (|2.ip lead to the 

following equivalent system: 

(2.3) 

r dtCi - Gi[Ci,c2](^i,V'2) = 0, 

dtC2 ~ G2[C2]V'2 = 0, 



dt(^P2dx1p2 - PlH[ClX2]{i^l,'ip2)j + g{P2 - Pl)dxC2 
+ ^5,(p2|9.V'2p-pi|i?[Cl,C2](^l,V'2)|- 



dxC2 



with 



P2(G2[C2]^2 + {dxC2){dxi^2)T " Pi(G2[C2]V'2 + (^.(2)^ [Cl , C2] (V'l , V'2)) ' 

2(1 + |5,C2P) 



2.2 Nondimensionalization of the system 

The next step consists in nondimensionalizing the system. The study of the linearized system 
(see 133], for example), which can be solved explicitly, leads to a well-adapted rescaling. Moreover, 
it is convenient to write the equations in the frame of reference of the ship. 

Let ai and a2 be the maximum amplitude of the deformation of, respectively, the surface 
and the interface. We denote by A a characteristic horizontal length, say the wavelength of the 
interface. Then the typical velocity of small propagating internal waves (or wave celerity) is given 

by ^ ^ ^ 

(P2 " Pi)did2 

' ; i — • 

P2di + pid2 

Consequently, we introduce the dimensionless variables 

^ _ Z ^ _ X + Cgt 

Z = — — , X = - , t 

di A 

the dimensionless unknowns 



Co = 



Co, 



C.(i) 



ijjiix) 



di 



ai ■ ■ a2Aco 

and the seven independent dimensionless parameters 
ai a2 



Pi 

P2 



£2 



di ' ^ di ' A^ ' (i2 ' Co' dicr 

With this rescaling, the full Eulcr system (12. ip becomes (we have withdrawn the tildes for the 
sake of readability) 

-eiFr9,Ci - -Gi (^1,^-2) = 0, 



di 

<5 = -i, Fr 

02 



Cs 

Co' 



Cn A^P2 
Bo = ^ ^ 



(2.4) 



{dt-FYd,K2 ~ -G2V2 = 0, 
P- 

(5t-Fr9,)(9,V2-7-ff(^i>2)) + (7 + '5)5.C2 
+ |9x(|9.V'2p-7|i^(V^i,V'2)|') - pt2d^M 



dxC2 



Bo nyi+T^^p:^;' 
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with 



(^G2^2+£2(9.C2)(g.V'2))' - 7(^G2V^2+£2(g.C2)g(^l,^2))' 
2(l+/i|e25,C2P) ' 

and the dimensionless Dirichlet-to-Neumann operators defined by 

Gi(V'i,V2) ^ Gf*[eiCi,e2C2](V'i,V'2) ^ -Aiei(^Ci)(5.0i) |..,+.,,, + |,^,^^^,^ , 



ax 

02^-2 = G^^'[e2C2]V'2 = -/^e2(a,C2)(a,(/.2) |_,,, + (9.</'2) | , 
H{ij,,ij2) = i?'^''[eiCi,e2C2](V'i,V'2) = +e2(a,C2)(a,0i)|_^,^ 

where and 02 are the solutions of the Laplace problems 

+ d^,) (1)2 = in n2 = {{x, z) e ]R2, -i < z < £2(2(2;)}, 

<i>2 = V'2 on r2 = {z = £2(2}, 

9^02 = onr6 = {z = -i}, 

+ dl) 01 =.0 in Oi = e R2,£2C2(a;) < z < l + £iCi(a;)}, 

01 = Vi on Ti = {z = 1 + £iCi}, 

9„2 0i = G2V'2 on r2. 

Remark 2.1. Let us keep in mind that in our case, the function Ci is not an unknown, but some 
fixed data of the problem: 

Ci{t,x) = Ci{x), 

where (^1 is the submerged part of the ship ( independent of time thanks to the change in the frame 
of reference). In that way, the first line of (|2.4I) is a relation connecting ipi with ip2, and (|2.4p 
can be reduced into a system of two equations with two unknowns. This reduction is computed 
explicitly in the following asymptotic models. 

Remark 2.2. The Cauchy problem associated to the full Euler system at the interface of two 
fluids is known to be ill-posed in Sobolev spaces in the absence of surface tension, as Kelvin- 
Helmholtz instabilities appear. However, in '341, Lannes proved that adding a small amount of 
surface tension guarantees the well-posedness of the Cauchy problem (with a fiat rigid lid), with 
a time of existence that may remain quite large even with a small surface tension coefficient, and 
thus is consistent with the observations. The idea behind this result is that the Kelvin- Helmholtz 
instabilities appear for high frequencies, where the regularization effect of the surface tension is 
relevant, even if Bo the Bond number measuring the ratio of gravity forces over capillary forces 
is very large^ On the contrary, the main profile of the wave that we want to capture is located at 
lower frequencies, and will therefore be unaffected by surface tension. Driven by this analysis, we 
decide to omit the surface tension term in the models of the Regime presented below, as they 
do not give rise to Kelvin- Helmholtz instabilities. On the contrary, we keep the surface tension 
term in the fully nonlinear models of Regime{^ As a matter of fact the regularization effect of 
the surface tension in this case plays a critical role in our numerical simulations, as it stabilizes 
the scheme, even z/Bo~^ is one order smaller than the coefficients in front of all other terms. 

The following terminology is used throughout the paper. 



^As an example, let us consider the values of the experiments of Vasseur, Mercier and Dauxois 1421 . One 
has dl = 0.05m, d2 = 0.12m, pi = 1 000.5kg.m~^, p2 = 1 022.7kg.m~'^. Taking advantage of the analysis of 
Lannes 1341 . we use a = 0.095N.m~^ as a typical value for the interfacial tension coefficient. As a conclusion, one 

has Bo = ""^^ — —. The coefficient in front of the surface tension term is therefore much smaller than the 
coefficients in front of any other term. 
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Definition 2.3 (Adapted solutions). We call adapted solution of (|2.4I) any strong solution 
(Ci: C2, V'lj "02), bounded in VF^'°°([0, T]; with T > 0, s > I and to big enough, and such 

that = Cii^) O'lT- the frame of reference of the ship) and the domains of the fluids remain 

strictly connected, i.e. there exists /imin > such that 

(2.5) 1 + eiCi (2;) - £2(2(^,2;) > /imin > and ^ + e2C2(i, a;) > Vin > 0. 

From the discussion above, it is legitimate to assume that such a smooth, uniformly bounded 
family of solutions of (|2.4p indeed exists. 



2.3 Description of the results, and the regimes under study 

The models displayed in the following sections are all justified with a consistency result; see 
definition below. When it is possible, that is when energy estimates are attainable, we provide a 
stronger justification with a convergence rate (the convergence theorems are precisely disclosed 
in Proposition 14. II and [ 



Definition 2.4 (Consistency). The full Euler system (j2.4p is consistent with a system of two 
equations (S) on [0, T] if any adapted solution U defines, via the changes of variables explained 
throughout the paper, a pair of functions satisfying (S) up to a small residual, called the precision 
of the model. The order of the precision will be 0(e), if there exists so,to > such that if 
U G W^''^i[0,T];H''+^°) with s > sq, then the L°°{[0,T]] H") norm of the residual is bounded by 
Co s, with Co a constant independent of e. 

Definition 2.5 (Convergence). The full Euler system (12. 4p and a well-posed system (S) of two 
equations are convergent at order 0{e) on [0, T] if any adapted solution with small initial data U 
defines, via the changes of variables explained throughout the paper, a pair of functions V such 
that V the solution of (E) with same initial data satisfies 



1^" ^L~([0,T];ff=) - '^0 ^' 



with Co independent of e. 



The small parameter e in these definitions is a function of some of the dimensionless parameters 
of the system that are assumed to be small. The regimes that we study throughout the paper 
have been briefly presented in the introduction; let us describe them here in more details. When 
nothing is specifically said in the regimes below, we assume the parameters to be fixed as 

7e(0,f), ei,e2e(0,f), < ^„,i„ < 5 < (5n.ax, Fr g [0, Fr„,ax], < Bo„,i„ < Bo . 
In particular, the two layers are assumed to be of finite, comparable depth. 

The mutual and crucial assumption between the two regimes is that the depth of the two 
layers of fiuids is small when compared with the internal wavelength. This is commonly called 
shallow water regime and is simply given by (with our notation) 

A* < f. 

The shallow water regime has been widely used in the framework of gravity waves. In the one 
layer, free surface case, it leads at order 0{ij) to the shallow water equations [31], and at order 
0{iJ?) to the Green-Naghdi equations ^5]. The analysis has been extended to the case of two 
layers with a free surface in [23] , and the models presented here can be recovered from models in 
this situation, when fixing the surface as some data of the problem. 
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Regime 1 (Small sized boats). 



M < 1 ; 



a EE ^ ^ ,1-7 



The two additional smallness assumptions of this regime are very natural. First, we assume that 
the depth of the submerged part of the ship is small when compared with the depth of the fluid, 
and the attainable size of the deformation. The numerical simulations of our model show that 
even with this assumption, the waves generated by the small disturbance can be very large (of the 
order of the depth of the two layers, and therefore very large when compared with the variation 
of the surface induced by the ship). This explains why the dead-water phenomenon can be so 
powerful. What is more, we assume that the densities of the two fluids are almost equal. This is 
known as the Boussinesq approximation, and is valid most of the time (the value of 1 — 7 reported 
in Celtic Sea in [SO] is ~ 10~^, and the value in the experiments of is w 10"^). 



Regime 2 (Small wave amplitudes). 



In this regime, we assume that the internal wave generated by the disturbance will remain small 
when compared with the depth of the two fluids. As previously, we assume that the waves 
generated by the ship are large when compared with the depth of the disturbance, so that the 
ship suffers from a relatively large wave resistance. In that way, the models we obtain involve 
only weak nonlinearities, and will be much easier to study. In particular, we are able to obtain 
well-posedness and convergence results for both of the models presented here. The counterpart is 
that these models remain valid only for relatively small waves. Again, this regime has been widely 
used in the literature. The one layer, free surface equivalent systems of the models presented here 
are the classical Boussinesq system [TUl E] and the KdV approximation [3D] . These models have 
been rigorously justified in [EH [71 [T], and the extension to the case of two layers, with a free 
surface, can be found in [35]. 



2.4 Expansion of the Dirichlet-Neumann operators 

The key assumption in our models, which is shared by both regimes [1] and [21 is the so-called 
shallow water assumption 



which states that the depth of the two layers of fluids are small when compared with the charac- 
teristic wavelength of the system. The main ingredient in the construction of asymptotic models 
in the shallow water regime, then lies in the construction and the justification of the following 
expansion of the operators Gi, G2 and H: 



Lemma 2.6. Let Ci, C2, V-i, ^-2, such that (Ci, C2, SxV'i, 9xV'2) e W^^°°{[G,T];H''+*«{W)), with 



o(m). 
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s > 1 and to > 9/2, and such that (j2.5p is satisfied. The 



has 



{T[h^,h2]d.^l +T[h2, 0]d,^2 



-i(/ii9|(/i25a;-02)) - {hiei-^Cidx{h2dxi^2)) 
2 dx 

G21P2 + ^.dx{h2dx''p2) - ti'^dxT[h2,0]dxi^: 



H 



(V'1,'02) - dxipi - fJ-dx(hi{dx{hidxi'i) + dx{h2dxi^2)) 

-\hldl^,^h,e,{dx^^)j^Ci{x)) 

with Co — Co ^7p— , l^lvi/i.oo/fs+to) ; '"^'^ operator T defined by 

T[h,b]V = ^^dxih^dxV) + ^{dx{h^{dxb)V)) - h\dxb){dxV)) + h{dxbfV. 

These estimates have been proved for L°°{[0,T]; H^) norms in Propositions 2.2, 2.5 and 2.7 
of [23]. The proof can easily be adapted to work with the time derivative of the functions, 
following the proof of [531 Proposition 2.12]. 



The idea is then simply to plug these expansions into the full Euler system (j2.4p . and drop 
all 0{^^) terms. Then, using the fact that Ci is a forced parameter of our problem, the system 
reduces to two evolution equations for (C2,f), with v the shear velocity defined by 



(2.6) 



V = ((02 - 70i) L-=,2C2 ) ^ 9a;?/'2 - 7-H"(V'l,1/'2)■ 



These calculations are postponed to Appendix |X] for the sake of readability, and we directly 
present here the system thus obtained: 



(2.7) 



idt-FTdxK2+dx 



hi + 7/12 



{hiv + jaFrC)) + fidx{C{vi,V2)) = 0, 



£2 , 



{dt-FTdx)v+{j + 5)dxC2 + ^dx 



\hiv + 7aFrCP - j\h2V - aFrCI^ 



+fJ-e2dx{Q[vi,V2]) = —d: 



{hi + 7/12) 



1 



dxC2 



v/1 + /iei|3xC2 



where C and Q are respectively linear and quadratic in (wi, V2), the latter being the approximation 
at order 0{fi) of {dxi'i,dxip2) given by 

h2V-aFT(i hiv + -faFTCi 

(2.8) vi = — ; — ; , and V2 ^ 



hi + 7/12 



hi + 7/12 



The operators £ and Q are precisely disclosed in Appendix [Xj The full Euler system (|2.4I) is 
consistent with this system at order 0{fi^) on [0, T], T > 0, as stated in Proposition lA.il page [331 

The simplified models wc study are deduced from system (|2.7p . using the additional assump- 
tions of Regime [Hand Regime [21 The full Euler system (|2.4I) is still consistent at order 0{ii^) with 
the models thus obtained, that is (see below) the strongly nonlinear systems (|3.ip and (|3.3p . the 
Boussinesq-type system (|4.2p and its symmetrized version (14.31) . The KdV approximation (|4.12p 
is then deduced from the symmetric Boussinesq-type system. 
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3 Strongly nonlinear models 

In this section, we introduce different strongly nonlinear asymptotic models for the full Euler 
system (|2.4p . These models are obtained starting from system (|2.7p . using the additional proper- 
ties of Regime \T\ The weakly nonlinear models (Regime [5]) are presented in following Section |31 
Indeed, the system (j2.7p only presumes some smallness on the shallowness parameter /i, that is 
to say the depth of the fluids are small when compared with the characteristic wavelength of the 
system. In the framework of Regime (TJ we assume that 

« 1 ; a = ^ = Oifi) , 1-7 = 0{p), 

where we recall that 7 is the density ratio and a the ratio between the amplitudes of the defor- 
mations at the surface and at the interface. As it has been said, page [Til these assumptions are 
very natural in the framework of our study. The first one supposes that the deformation induced 
by the presence of the ship at the surface is small when compared with the depth of the two 
layers of fluid, and small when compared with the attainable size of the wave generated at the 
interface. The second one is the classical Boussinesq approximation. The simplified systems we 
obtain remain consistent at order C'(^^), as stated in Proposition 13. II In Section [321 we present 
and discuss numerical computations from the constructed models, that reproduce the key aspects 
of the dead-water phenomenon. 



3.1 The fully nonlinear model in Regime [T] 

The first obvious observation is that the total depth of the fluid is approximatively constant in 
Regime [11 i.e. equal to h = 1 -|- | at order C'(//). What is more, one has 



1 + ^ + o{^l). 



hi + 7/12 ^ h + 0{fi) 
Therefore, the approximations of {dxipi-, dx'ip2) at order 0(/i) given in (|2.8p are now simply 

/12W — —h vi + 0{p), hiv = h V2 + 0{^i). 
It follows then some substantial simplifications, and one obtains in the end the system 



(3.1) 



{dt^Fvdx)C2 + 



hih2 



a Fr 



hi + 7/12 



Ci) + t^dx{Viv) = 0, 



(-i V \ , / , A , ?i ,'1 Ift-ii'P - 7|/i2Wp aFr 



2 {hi+^h2f 



1 



l^t2dx{V2[v]) = —dl 



dxC,2 



where h = 1 + h and the operators Vi and V2 are defined by 



h'^VlV = -{h2dx{hldx{h2v)) + hidx{hldx{hiv))), 
h^ V2[v] = [hldx{h2v) - hldxihiv)] + i {{hidx{h2v)f - {h2dx{hiv)f) 
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Linearizing the system p.ip around the rest state, leads to 

1 „ 1 



{dt^FYd^)v 



7 + 5 



dxV 



35(1 + 5) 



a Ft 



dxCi, 



from which we can easily deduce dispersion relation of the system, without forcing. Indeed, setting 
a = 0, the wave frequency a;(fc), corresponding to plane- wave solutions g^^-^-*"*^ jg solution of 
the quadratic equation 



(3.2) 



(Frfc + ojY = { 1 



Bo(7 + 5) 



2_S+J_ 

3(5(1 + (5) 



Therefore, for high wave numbers, k, there is no real- valued solution Lu{k) of (|3.2I) . which means 
that the system p.ip is linearly ill-posed, and leads to instabilities. In order to deal with this 
issue, we use a nonlinear change of variables on the shear velocity u, that leads to a formally 
equivalent system. This idea is not new (see [51 SHI HSl US] for example), and usually relies on a 
"shear velocity" constructed from horizontal velocities at any depth in the upper and lower layers. 
The choice of the depth allows then to control properties of the linear relation dispersion. Here, 
we define w with 



hih2 



hih2 



hViw 



/1I+7/12 hl+^h2' ■ ^~ ■ ~ r- (-/j_^2)/l2 

This leads immediately to the following system, equivalent at order 0{fi^): 

hih2 



(3.3) 



idt-FTdx)C2 + dx 



h2^CT\ =0, 



hi + 7/12 h 
{dt ~ Fr dx) (w - ^iSiw) + {-/ + S)dxC2 + ^28, 



fie2dx {w S2W) = -^dl 
Bo 



1 hj-jhl 

2 (/ii +7/12)2 

1 



a Ft 



dxC,2 



\P^+ij^\dIC2 



with /i = \ + \,hi = 1 + eiCi - e2C2, ^2 = 



e2C2, and the operators 



SiW= -{dl{hih2w) ^{dxh2Yw) 



S2W = 



hih2 
~ih 



{{dlh2)w + 2{dxW){dxh2)) 



hi - h2 
2h 



{dxh2fw. 



Under this form, our system corresponds to the Green-Naghdi model presented in the one-layer 
case in [3 [5], and proved to be well-posed and convergent, in the sense that their solutions provide 
an approximation at order ©(/i^) in L°°([0, T]; H^) to the solutions of the one layer water wave 
equations. When dropping all 0(/i) terms, one obtains a shallow water model, similar to the 
ones derived in [121 and ^20^ in the flat rigid lid case (the forcing terms induced by the presence 
of the body do not appear). Such a system has been studied in details and analyzed as a 
hyperbolic system (leading to well-posedness results under reasonable assumptions on the initial 
data) in [H [271 [13 . 
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One can check that the issue of hnear ill-posedness is now solved. Indeed, hnearizing the 
system (|3.3p around the rest state, leads to 

1 aFr 

{dt-Frdx)C2 H r^d^^w = — — r^^^Ci, 

7 + d 1 + 

(5,-Fra,)(l-|^a> + (7 + 5^2 = ^9^C2. 
The dispersion relation of this linear system, when setting a = 0, is 

which leads to real- valued waves frequencies uj, for any values of fc G M. 

Let us remark that under the assumptions of Regime [U the relations p.2|) and (j3.4p are 
asymptotically equivalent at order ©(/x^A:^), so that the effect of the nonlinear change of variable 
only affects high frequencies. 



We now state that the two systems (j3.ip and p.3p are equivalently justified as models for the 
full Euler system, with the following consistency result. 

Proposition 3.1. Assuming that a = 0{fi) and 1 — 7 = 0{fj,), the full Euler system (|2.4p is 
consistent with the models (13. ip and (j3.3p . both at precision 0{^'^) on [0,r], tuith T > 0. 

Proof. Let U = (Ci, C2, "01, -02) be a strong solution of ([231), bounded in W'^'°^{[0, T]; with 
s > 1 and to > 9/2, and such that (|2.5I) is satisfied with (i{t,x) = Ci{x). The consistency 
result of Proposition lA.ll states that {(2,^), with v = dxip2 — (V'l; 02), satisfies (|2.7p , up to 

Ri = iri,r2)^ e L°°{[0,T];H')^, with (for i = 1,2) 

2/1 

Since TZi, TZ2, 7" and dxT-L (defined in Appendix VK\ involve two spatial derivatives, and thanks 
to straightforward calculations using the smallness assumptions of Regime [TJ one has for i = 1, 2 



\L, - V, 



L°=([0 



\i ''mm / \ ' 'mm / 



The consistency of (|2.4p with p.ip is therefore proved. 

Now, we set w = v + M Pi [/t2]f , and one has immediately, using the fact that H^iM.) 

is an algebra for s > 1/2, 



hVi [h2]'w 

V — w + fi 



{h - h2)h2 



< fi^ Co 



- '•'^ ( ir~' l^lwi ^ff^+'o 



It is then straightforward to deduce from the previous consistency result, the consistency of (|2.4p 
with ([33]). □ 
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3.2 Numerical simulations 

In figures [3] and m we plot the behavior of the flow predicted by the model (|3.3p . for different 
values of the parameters. Each figure contains three panels. The left panel represents the interface 
deformation, depending on space {x € [—20, 20]) and time {t e [0, 15]) variables. The right panel 
contains the time evolution of the wave resistance coefficient, computed thanks to formula (|C.4|) . 
Finally, we plot in the bottom panel the situation of the system (i.e. the surface and interface 
deformations) at final time t — 15. 

In these simulations, and throughout the paper, we use zero initial data conditions, and a 
function Ci{x) defined by 



The parameters of the system 62, a = £1/^2, S, 7, Fr and Bo, are specified below each of the 
figures. The schemes used are described and justified in Appendix ID] 

Here, we decide to study the effect of the depth ratio coefficient d. We choose two different 
values for the depth ratio: i5 = 5/12 (corresponding to the experiments of [H]), or i5 = 12/5 for 
a thicker upper layer. We plot in figures [3] and |3] the outcome of our scheme (|D.4I) . setting the 
Froude number as Fr = 1. Away from this critical value, as predicted by experiments |25ll42j and 
confirmed by our numerical simulations, the amplitude of the generated waves (and therefore the 
magnitude of the wave resistance coefficient) are significantly smaller. In that way, the outcome is 
similar to the one predicted by the weakly nonlinear models, and we do not present these results 
for the sake of brevity. We let the reader refer to Section 14.3.11 for a discussion on the case of 
subcritical and supercritical Froude values, in Regime [H 

Let us express a few remarks on these results. There are some similarities for each of the 
situations. First, one sees that there is no deformation on the right-hand side of the plots, which 
corresponds to the up-stream propagating part. This is due to the fact that we have set Fr = 1, 
which corresponds to the maximum gravitational wave velocity in the flat rigid lid case. On the 
contrary, on the left-hand side (down-stream propagating part), one remarks a small elevation 
wave, followed by even smaller perturbations, progressing with velocity c_ = —2 (in the frame 
of the ship). This corresponds to the rj- part of the KdV approximation decomposition, and is 
studied with more details in Section H31 Finally, the most important part takes place just behind 
the location of the body. Each time, an important wave of elevation is generated just behind the 
body, producing a severe wave resistance (see Remark FC . 1 1 page HO]) . This wave comes with a tail 
of smaller waves, that are located away from the boat, and therefore do not produce any drag. 

However, one can see that the shape and time-behavior of this elevation wave is quite different, 
depending on the value of the depth ratio S. When the upper fluid domain is thinner than the 
lower fluid's, the generated wave is flattening as it is growing up. Its height is relatively low, 
but the deformation carries on with a depression wave, located just below the body. On the 
contrary, when the upper fluid domain is the thicker one, then the generated wave remains sharp 
and is continuously growing. It is separated with its tail, while all of the energy produced by 
the ship contributes to the elevation of the wave. The wave produced wave resistance is greater 
in this case, and the dead-water effect appears therefore much stronger. The fact that the wave 
resistance is stronger when the upper fluid domain is thicker is counter-intuitive, as the distance 
between the ship and the interface is larger. However, as we can see in Appendix ([Cl, the wave 
resistance coefficient does not depend on the distance between the two interfaces (Ci — ^2), but 
rather on the amplitude of the generated wave. It is therefore not surprising that, when the 
upper layer is thicker, the generated elevation wave is able to reach larger amplitudes, and the 
drag suffered by the body is stronger. 




ifxe (-1,1) 



otherwise. 
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Figure 4: Flow predicted by (I3.3|) . with steady initial data and a thicker upper layer. 
€2 = 1, a = Ai = 0.1, 7 = 0.99, Fr = 1, Bo = 100, 6 = 12/5. 
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Figure 5: Flow predicted by p.3p . when the velocity of the body, Fr, is not constant. 
€2^1, S ^1, a = fi = 0.1, 7 ^ 0.99, Bo = 100. 

The experiments conducted in [55] (this has also been reported in HU] for example) exhibit 
an interesting behavior, that cannot be seen in our simulations in figures [3] and 2) Indeed, with 
their setting, the dead-water phenomenon appears to be periodic in some sense, as the generated 
wave increases its amplitude, slows down the boat, breaks, and the process repeats. As a matter 
of fact, even when running simulations for much longer (around ten times) time intervals, the 
solution of p.3p does not display such a phenomenon. It is interesting to see that this is a major 
difference with the solutions of the KdV approximation, which generate periodically up-stream 
propagating waves, inducing an oscillation in the related wave resistance (see Figure [91 page [50)1 . 
However, even with this model, the time-period of generation of these solitons (AT ss 100 in our 
situation) is too large to explain positively the phenomenon. 

This discrepancy may be due to our assumption of a constant velocity for the body, as the 
setting in [42) constrains a constant force to move the boat. As a numerical experiment, we 
performed simulations with a Froude number Fr adjusted at each step: 

Fr((n + l)At) = Fr(nAt) - AtCstti{Cw{nAt) - Cstt2). 

This roughly corresponds to a case where the acceleration of the ship is given by a constant 
force imposed to the ship, minus the resistance suffered by the bodyn We present in Figure [5] 
the result of this computation. A periodicity can clearly be seen in the wave resistance (with a 
time period of order AT « 10). One can explain the phenomenon as follows. The velocity of 
the ship, suffering from the drag generated by the wave resistance, decreases down to a value 
where the generated wave resistance is very small. Released from its drag, the ship speeds up, 
and the phenomenon repeats periodically. This corresponds to the observations of [42], and is 
not recovered by constant velocity models (see [14, 45, 59 ). 

^Of course, this modification is careless, and the produced scheme reflects only roughly what would be a model 
where a constant force is imposed to the body. To construct and rigorously justify a constant-force model would 
be much more difficult, as we explicitly use the constant- velocity hypothesis, as early as in I IA.41 1. 
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4 Weakly nonlinear models 

The aim of this section is to introduce simple weakly nonlinear models for our problem, using 
the smallness assumptions of Regime [5J 

« 1 ; 62 - 0{fi), a = 0{fi). 

We recall that /i stands for the shallowness parameter, that is the ratio between the depth of 
the upper fluid layer and the internal wavelength. The parameter £2 measures the magnitude of 
the deformation at the interface when compared with the upper fluid layer, and a is the ratio 
between the amplitudes of the deformations at the surface and at the interface. As discussed 
in Remark 12.21 the effect of the surface tension term will be relevant only if Kelvin-Hclmholtz 
instabilities, located at high frequencies, appear. The models we obtain in this section do not give 
rise to these instabilities, and we therefore neglect the surface tension term, and set Bo^^ = 
(Bo~^ = 0{fi'^) would in fact suffice). 

Starting with the strongly nonlinear model (12. 7p . we deduce easily a Boussinesq-type model 
and its symmetrized version. Then, using a classical WKB expansion, we obtain a rougher 
approximation, that consists in two uncoupled Korteweg-de Vries equation, with a forcing term. 
These two equations are studied in details, and numerically computed, in Section 14.31 The 
symmetric Boussinesq-type model, as well as the KdV approximation, are justified thanks to 
convergence results (see Propositions 14. II and l4.2l respectively). 



4.1 The Boussinesq-type models 

Let us drop 0{fP) terms in system (j2.7p . using the assumptions of Regime [5] One obtains 

straightforwardly the following system 

(4.1) 

idt-Frd,)v + h + S)d^C2 + '-l-f—2_d4\v\^) = 0. 

This Boussinesq-type system can be written in a compact form as 

(4.2) dtU + Aod^U + €2Ai{U)d^U - ^iA2^lU = a&o(a;), 

withC/ - (C2,«)^, &o = -Fr3^(£Ci,0)^, and 

Following the classical theory of hyperbolic quasilinear equations, our aim is now to obtain an 
appropriate symmetrizer of this system. Let us define 

S{U) = 5o + e2Si{U) - fiS2dl 

Multiplying (j4.2p on the left by S{U) — ^KSod^, and dropping the 0{fi^) terms, leads to the 
following equivalent system 



(4.3) (^0 + e2Si{U) - fi{S2 + KSo)dl)dtU -f ( So + e2^i{U) - fi^2 ) d^U = ab{x), 
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with 6 = -Fr7(£Ci,0)^, and 



i + is fo i\ , ^ 



35(^ + 7) VI 

As we can see, the system (I4.3P is perfectly symmetric, and K can be chosen so that S'o and 
S2 + KSq are definite positive. With these properties, we are able to use energy methods in order 
to prove that (j4.3p is well-posed, and convergent with the full Euler system (j2.4p . at order 0{iJ,) 
up to times of order 0{l/ii). More precisely, one has 

Proposition 4.1. Let s > 3/2, to > 9/2 and U — (Cii 5 "015 "02) be an adapted solution of the 
full Euler system 1^^, bounded in M/i^°°([0, T/^); We define V = (C2,w) by 

V= 5:r (02 - 701 ) = 9xll^2-jH{lpi,1p2)- 

Moreover, let us assume that there exists a constant Cq such that 62 < CofJ- and a < Cofi. Then 
there exists a constant Ci — C{::^,j + 6,Co) > such that if e2\V \t=o — c^T' ^^^'^^ 

ists T > 0, independent of ^, and a unique solution Vb G C^IP, T//i); n C^([0, T//i); if^), 

hounded in W^^'°°([0, T/^); iJ^), of the Cauchy problem with Vb Uo = V . Moreover, 

one has for all t G [0,T/fi), 

Before starting with the proof, let us remark that the proposition is not empty only if there 
actually exists a family of solutions of (I2.4p . smooth and bounded in W^'°°{[0,T/ fi); H^~^*°). As 
discussed in Remark l2.2l page[^ this requires adding a surface tension term. However, this surface 
tension term is very small in practical cases, so that we can assume ^ = ©(/i^), and the result 
is obtained as in the proof presented below. 

Proof. Step 1: Well-posedness. In order to prove the well-posedness of the symmetric sys- 
tem (|4.3p . we use techniques of It is proved in Proposition 2.4 therein that systems of 
the form (|4.3D (with four equations instead of two and without the right-hand size) , satisfying 

i. The matrices S'o, Soi '5'2, '^2 are symmetric, 

ii. Si{-) and Si(-) are linear mappings, and for all U E M*, Si{U) and are symmetric, 

iii. S'o and S2 are definite positive, 

are well-posed and satisfy an energy estimate. The proof is easily adapted for the case of a non- 
zero right-hand side, and we briefly give the arguments here. The key point of the proof relies 
on a differential inequality, satisfied by the following energy (with A* = (1 — d'^)'^^^) 

Es{U) = l/2{SoA'U,A'U) + e2/2{Si{U)A'U,A'U) + ti/2{S2A'd^U,A'd,U). 

This energy is easily proved, thanks to the positiveness of Sq and S2 and using the smallness 
assumptions of Regime[51 to be equivalent to the I • I 3+1 norm, that is to say there exists Co > 
such that 
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Then, we prove that there exists Ci — C'i(:^^,7 + i^) > such that if e2|C^|^s+i < then the 

operator defined by P{U,dx) = So + e2Si{U) — f^S2dl : H"^^ — >• if""^ is one-to-one and onto, 
and that P{U,dx)~^{^o + e2^i{U) — ^S29^) is uniformly bounded — H^, so that any 
solution V of (14. 3p will satisfy the a priori estimate 



\dtV\ 



H; 



with C2 independent of £2 and /i, as long as e2|^|^s+i < It follows that Es{V) satisfies 

^Es{V) = e2/2iSiidtV)A'V,A'V) - e2{[A' , Si{V)]dtV,A'V) + e2/2iiJ:,{d:,V)A'V),A'V) 
at 

-e2([A^Sl(F)]a,F,A^y) + a(A^6,A^\^) 
< C4 [e2Es{Vf^^ + a\b\^^Es{Vy 



Nl/2 



The last inequalities come from Cauchy-Schwarz inequality, Sobolev embeddings, Kato-Ponce 
commutator estimates and the above a priori estimate (see Appendix A of [22 ). The Gronwall- 
Bihari's Lemma allows to conclude that as long as e2|V|^a+i < one has 

(4.4) < Co E{Rsf'^ < C5a|/|^, t. 

Using the assumptions of Regime [2] (a — 0{fi),e2 — 0{fi)), one can then follow the classical 
Friedrichs proof, and obtain the existence of r(7 -I- S, {^2\V^\^s+i)^^) > and a solution V 

of g31), defined over times [0,T/^), such that V € C°i[0,T/ fi); H'+^) n C^{[0,T/ fi); H"). 

Finally, the uniqueness of the solution is obtained in the same way, applying the energy esti- 
mate to the difference of two solutions. Indeed, let ^1,^2^ C°([0, T/e); H^+^) n C^{[0, T/e);H') 
be two solutions of the Cauchy problem (|4.3p with same initial value Vi \^^g — V2 \^^g — V'^. The 
functions Vi and V2 are uniformly bounded thanks to (|4.4p . One can immediately check that 
R=Vi—V2 satisfies 

(^SQ + e2S,iVi)-fiS2dl)dtA'R+ (^^0 + ^2^iiVi) - ^i^2^l)^^A' R 

(4.5) +e2[A%Si{Vi)]dtR + e2[A^j:,{Vi)]dxR^e2F, 

with F = — A*^S'i(i?)i9(V2 -I- 12i{R)dxU2^ ■ Then, we can carry out the above method on R with 
a modified right-hand side, and obtain the equivalent energy estimate 

J^Es{R) < e2CG{\Vi\jj^ + \V2\j^^,)E,, 

with Cg = C'(:^^, i5 + 7, l^'^l^f+O- From Gronwall-Bihari's inequality, using the uniform bound- 
edness of Vi and V2 on [0,T/e), and since Es{R) I^^q = 0, one has immediately Es{R) = on 
[0,r/e), and finally Vi = V2. 

Step 2: Consistency. We prove that, assuming that a = ©(/i) and €2 = 0{ji), the full Euler 
system (|2.4p is consistent with the models (j4.2D and (j4.3p . both at precision 0{fi^) on [0,r], 
T > 0. 

Let U EE (Ci,C2,V'i,V'2) be a strong solution of (EH), bounded in W^^°°{[0,T]; H'+*o) with 
s > 1 and to > 9/2, and such that (|2.5p is satisfied with (i{t,x) = Cii^)- The consistency 
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result of Proposition lA.ll states that (C2,'^), with v = dx^J2 — 7^^(V'ii V'2), satisfies (|2.7p . up to 
Ri = iri,r2f e L°°{[0,T];H')^, satisfying (for i = 1,2) 



1 



It is the straightforward, as in the proof of Proposition 13.11 to show that the terms dropped 
in gH, and later in (g^]), are ah bounded by /i^ Co (7^, |C/|^i ,^^,+4) in L~([0, T]; 7?") 
norm. Finally, {C2,v) satisfies (|4.2I) and (14. 3p with modified right-hand sides, i.e. , respectively, 

abo — abf) + /i^/o, and ab = ah + ii^f, 

with /o, / uniformly bounded in i°°([0, T];H'). 

Step 3: Convergence. The convergence estimates of the proposition follow easily from the cal- 
culations of Step 1, together with Step 2. Indeed, thanks to the consistency result, V — Vb 
satisfy (|4.5p . with modified right-hand side 

F = F + fi^f, f eL°-{[0,T];H'). 

It follows that Es{V — Vb) satisfies 

j^EsiV-VB) < e2CeE,{V - Vb) + fi^Cjlfl^^EsiV -Vb^^^ 
and the Gronwall-Bihari's Lemma leads to 

\y-yB\,^ao^T/,y,H^.^^) < CoE,,iV-VB)'^' < g|/|^.(e^-*-l). 
The convergence estimate of the proposition follow from e2—0{p). □ 
4.2 The Korteweg-de Vries approximation 

In this section, we use a WKB (Wentzel-Kramers-Brillouin) approximation, in order to deduce 
from the symmetric Boussinesq-type system (j4.3p . an approximated model that consists in two 
uncoupled forced Korteweg-de Vries (fKdV) equations. This method has been used, for example, 
in [22I [35] , and is briefly presented in the following. 

The idea is to look for an approximate solution of the Cauchy problem (j4.3p with initial data 
C/^, under the form 

C^app(i,a;) = Uo{^it,t,x) + ^iUi{fit,t,x), 

with the profiles Uo{T,t,x) and fiUi{T,t,x) satisfying Uq It^^^o = C/° and Ui l^^^^o — 0- 
Plugging the Ansatz into (14. 3p leads to the following equation 

{Sodt + ^od.)Uo + fiSodrUo + e2{Si{Uo)dtUo + Si(C/o)a,C/o) - ^i{S2dldtUo + Esa^f/o) 

(4.6) +/i(5o9t + Soa,)C/i + Ai'i? = 0. 

We now deduce Uo{T,t,x) and Ui{T,t,x), by solving (|4.6p at each order. 
At order 0(1) : We solve 

(4.7) ( Sodt + )Uo = 0. 
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Let US define e± = (zfc , Vl + ^)'^ ■ One can eheck that the basis satisfies 

ei ■ SoBj = CiSij, and • SoBj = 5ij, 

with 6ij the classical Kronecker delta symbol, and c± — ±1 — Fr. 

Therefore, when we define u± = e± ■ SqUq (and hence Uq — u+e+ + ii_e_), the scalar 
product of (|4.7|) with e± leads to 

( dt + c±dx )u± 0. 
Finally, since u± satisfies a scalar transport equation, we use the notation 

(4.8) u±{T,t,x) — u±{T,x — c±t) = u±{t,x±), 
with initial data u±(0, x±) e± ■ SqU'^{x±). 

At order C'(^) : We solve 

(4.9) S„drU„ + - (Si([/o)a,C/o + Si{Uo)dtUo)-^2dlUo-S2djdtUo + {Sodt + J:od.)Ui = -b{x), 
that we can split in 

(4.10) drU± + X±u±dx±u± + iy±dl^u± = l3±(x), 

with A± = — e± • (Si — c±S'i)(e±))e±, = e± • (— E2 + c±S'2)e± and /3± = — e± • 6; and on 
the other hand, 

(4.11) {dt + Cidx)e.r SqUi + ^ aijkUk{T,x - Ckt)dxUj{T,x - Cjt) l3ijdluj{T,x - cjt), 

with a^jk = ei ■ (Si(efc) - CjSi{ek))e-i and /3y = • (S2 - CjS2)ej. 

It is clear that Ui satisfies (|4.8p and (|4.10p . if and only if Ui{et,t,x) satisfies the Korteweg-de 
Vries equation: 

dtu± + c±dxU± + /i ( A±M±(?2;U± + v±dlu± ) = M/^±(a;)- 
Finally, simple calculations show that in our case, we can decompose 

C2 = ?y+ + with ?7± = ± , / =M± 

^2(7 + 5) 

satisfying precisely the following KdV equation 

dtv^ + (-Fr±l)a.,;± ± ,2^^r;±9.^± ± = -"Fr^^Ci- 

Unsurprisingly, we recover the KdV approximation with a flat rigid lid, when a = (see [22j 
and references therein). 



Using these forced Korteweg-de Vries equations as an approximation of the full problem is 
justified up to times of order 0(l//i) by the following proposition. 
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Proposition 4.2. Let s > 1/2, to > 5 + 5/2 and U = (Ci; C2, V'l: V'2) &e an adapted solution of 
the full Euler system (HH), bounded m VKi'°°([0, T/^); We define V = (C2,v) by 

V = (02 |, = ,2C2 ~ 701 L- = ,2<:2 ) = 9a;V'2 - T-ffCV'l, V'2)- 

Then there exists 77+ and the two solutions of the following forced Korteweg-de Vries equation 

(4.12) dtv± + (-Fr±l)a,7;± ± e2~^V±d^ri± ± ^x^j±^dlrj± = -aFT-f-^C,{x), 

with rj± Ij^j, = \ {C,2 ± ^^^^''^) lt=o • Moreover, if there exists a constant Co such that £2 < Co/i and 
a < Co/x, i/ien one ftas /or all t € [0, T/^), 

IC2 -(??++ ?7-)|i„([o ,].j:,, + |«-(7 + '5)(77+-?7-)Lo„([o^,].^.) < A'VtC, 

with C = C(^, 7 + S, |C/|^i,oo^,+to ' C'o)- 

Additionally, ifV satisfies (1 + a;^)V"|j^„ G iJ""^^, i/ien one /las t/ie improved estimate 

|C2-(?7++?7-)L^([o^,].^e + |^^-(7 + '5)('/+-'7-)L»([o,t];^^.) < eC^', 

C = C(^, ^,7 + -5, |f/|^yi.„„^,+,o , 1(1 + x^)V 1^,+,). 

The proposition is obtained by a simple adaptation of the techniques presented in [22], with 
additional forcing terms in the Korteweg-de Vries equations. The proof is given in Appendix [Bl 
for the sake of completeness. 

4.3 Analysis of the forced Korteweg-de Vries equation 

The forced Korteweg-de Vries equation 

(4.13) dtu + cdxU + eXudxU + evd'tu = e-^f(x) 

ax 

has attracted a lot of interests, especially in the framework of the one layer water wave problem 
(where a moving topography, or pressure, is the forcing term that generates waves). Of particular 
interest is the problem of the generation of solitons, that have first been numerically discovered by 
Wu and Wu , and validated with experiments by Lee [35] . Using the Boussinesq-type system 
or the KdV approximation, they found that starting with a zero initial data, the solution can 
generate periodically an infinite number of solitary waves. Numerous work have have then tackled 
this problem, including [SHI ISZl liSl ISI] ■ It appeared that the Froude number (which is given by 
1 — c in (j4.13l) ) is playing a predominant role in this phenomenon, as the generation of solitons 
only occurs for a narrow band of its values. One could roughly summarize the observations by 
the existence of Fc > 1 such that 

i. if Fr > Fc, then the flow approaches a steady state, symmetric and localized at the site of 
forcing; 

ii. for Fr < Fc, solitons are periodically generated at the site of forcing and radiated up-stream; 

iii. the amplitude of the generated solitons goes to zero as Fr — > ~oo. 
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The existence of steady solitary waves of (|4.13p . and their stabihty is therefore essential. This 
issue has been studied for specific forcing terms by Camassa and Wu in [T3J [T3] , and for general 
forcing terms by Choi, Lin, Sun and Whang in [T71[TH]. They prove that for Fr sufficiently large, 
there exists a unique small steady solution. Moreover, this solution is proved to be symmetric 
and localized at the site of forcing, and stable in iJ^(R) (in the Lyapunov sense). The behavior 
when the Froude number Fr approaches unity is more peculiar and depends on the sign of the 
forcing term. Subcritical Froude numbers are not studied. 

In the framework of our analysis, the values of the coefficients A and v depend on the param- 
eters 7 and (5, and their order of magnitude depends on ti and /i (the magnitude of the forcing 
term depends mostly on a). We use the factor e to keep in mind that all of these coefficients are 
assumed to be small. More precisely, we recall 

3 J2 - 7 1 1 + 7^ 

c± = ±1-Fr, eA± = ±£2- ^ _^ j : e^^i = =^^6 5(7+7)' ^^'^^^ " - Fr 07(1 (x) . 

When (5^ — 7 ^ 0, a simple change of parameters allows to recover values of the one layer problem, 
and one obtains straightforwardly similar results. Our approach is quite different. In the following 
section, we prove that away from the critical speed (c = 0), the solution of (|4.13p with small initial 
data will remain small (in iJ*(R) norm) for long times, using smallness assumptions on the forcing 
terms and the coefficients eA and ev. Numerically, it appears that the solution converges locally 
towards a negative steady state in the supercritical case, and that small solitons are continuously 
generated otherwise. Then, in Section 14.3.21 we study numerically the generation of up-stream 
propagating solitons for Froude numbers around the critical value. 

The consequences of this study to the dead-water effect are the following: 

i. away from the critical Froude number, the drag suffered by the ship is always small; 

ii. the peak of wave resistance occurs for Froude numbers just below the critical value; 

iii. the time-period of the generation of solitons predicted by the KdV approximation is of the 
same order as the time-scale of relevance to the model, so that in cannot clearly be held 
responsible for the periodic aspect described among others in [12] . 



4.3.1 Non-critical Ftoude numbers 

In the following proposition, we obtain an improved growth rate for the solution of the forced 
Korteweg-de Vries equation (j4.13p . if the velocity coefficient c is away from zero, when assuming 
smallness on the nonlinearity and dispersion coefficients A and v, and on the initial data and the 
forcing term. This phenomenon is easily explained, when looking at the linear transport equation 
related to (|IT^ : 

(4.14) dtv cd^v = £^/(2:). 

The Cauchy problem, with initial data v\^^^ = euq, is solved by 

V = e uo{x - ct) + ^ f{x) - f{x - ct)y 

The solution is therefore bounded for all times as soon as c 7^ 0, and small if the initial data and 
the forcing terms are small. We will obtain improved bound estimates on u the solution of the 
forced KdV equation, by estimating the difference kt — w 
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Proposition 4.3. Let s > 3/2 and u be the solution of 

dtu + cdxU + eXudxU + evd^u — e—f{x), 

dx 

such that M 1(^0 = suq. Let us assume that there exists M > such that 



Then there exists T{M,s) > and C = C(M, s) > 0, such that the function u is hounded on 
[0,r/£] by 

l"WL~([0,T/e];//=) - 

Proof. Thanks to [5], we know that the function u exists and is unique. We define as above 

V = e uoix -ct) + ^ f{x) - f{x - ct)^, 
the solution of the transport equation 

dtv + cdxV = ^^■^f^^)' 
with same initial data w Ij^q — suq. Since the function v is uniformly bounded by 

the estimate of the proposition follows from the same estimate on the difference r = u — v. It 
is easy to check that r satisfies 

(4.15) dtr + cd^r + eX{r + v)dx{r + v) + evdl{r + v) ^ 0, 

and r \^^^ = 0. When multiplying (|4.15p by A^'^r, and integrating, one obtains 

i^(AV,AV) + c(A^a,r,A^r) + e\{k' {{r + v)dx{r + v)) , r) - ey{k'dl{r + v),K'dxr) = 0, 

denoting (/, g) = Jj^ fg the L^-based inner product. It follows then 

\j^{V\\) = -e\{K^{r + v)dx{r + v)),K'r) - ei^{K' dlv, K'r). 

Some of the terms of the right-hand side are straightforwardly estimated using Cauchy-Schwarz 
inequality, and the algebraic properties of iJ*(R), s > 1/2: 



V 
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2 


V 




r 





As for the remaining terms, we integrate by part and obtain 

(A^(ra,r), AV) = -i(9,rAV, AV) + ([A^ r]9,r, AV), 
(A^(z;a,r),AV) = -i(9,«AV, AV) + ([A^ i;]a,r, AV), 
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with [T, f] the commutator defined by [T, f]g = T{fg) - fT{g). Since for f £ H', s> 3/2, one 
has |c)a;/|^^ < one can use the classical Kato- Ponce Lemma [55] to deduce the following 

estimates 

\{A^vd,r),A'r)\ < Cf,{s)\v\^^\r\l^, 
|(A^(r5,r),AV)| < Co{s)\r\l^. 

Altogether, it follows that \r\^^ satisfies the differential inequality 

1 d /| |2 \ ^ ^ , , , /| |2 II II I |2 I |3 \ II II 

with Co(s) a constant depending only on the parameter s > 3/2. Now, using the assumptions of 
the proposition, it is straightforward to see that one can rewrite 

Now, let us define T* the maximum time such that \r\^^, < e on [0,T*] (this is true at i = 0, 
so that r* > by a continuity argument). Then we deduce from the above inequality, and 
Gronwall's Lemma, that for all t £ [0,T*], 

with a constant C(s) depending only on the parameter s > 3/2. From this very a priori estimate, 
we deduce that T* > (C(s)e)~^, and the proposition is proved. □ 

Remark 4.4. One of the immediate consequences of Proposition \J7S\ to our problem is that, in the 
decomposition of the KdV approximation ( Proposition \4-'^ , the function rj^ remains small on the 
time interval where the KdV approximation is a relevant model. Indeed, as c_ = —1 — Fr < —1, 
then if rj^ = (small in H"^^^ would suffice), then for t £ [0,T/fi), one has 

Therefore, we focus in the following on 77+, the solution of 

(4.16) dtv+ + il-Fr)d,V+ + ^4^-^^+^-^+ + ^'lJ(^^)^''^+ = -"^'^i^i- 

In Figures |6] and [71 we compute the solution of (j4.16p . with zero initial data, in the supercritical 
(Fr = 1.5) and subcritical (Fr = 0.5) cases. For each of the figures, we plot the flow, depending 
on space and time variables, as well as the evolution of the related wave resistance coefficient Cwi 
calculated with formula (|C.5|) page[ini We use two values for the depth ratio: 6 £ {5/12, 12/5}. 

As we are away from the critical Froude number, the amplitudes of the generated waves, as 
well as the magnitude of the wave resistance coefficient, are small (at most of the order 0{fj,)). 
One remarks that the behavior of the flow is roughly independent of the value of the depth ratio 
6. This is easily explained, as the quadratic nonlinearities do not play an important role when the 
deformation is small (as they are of order 0{fj,^)). Therefore, the variations of 5 are mostly seen 
through a variation of the size of the dispersion coefficient v, which do not change the behavior 
pattern of the solution. 

On the contrary, the Froude number parameter has an essential role on the behavior of the 
solution. One sees that for subcritical fiows (Figure [5]) , up-stream propagating solitary waves are 
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Figure 7: Supercritical flow: Fr = 1.5. a = e2 = /x = 0.1, 7 = 0.99. 
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continuously generated, as a widening oscillatory tail is generated at the down-stream propagating 
area. The former will generate most of the wave resistance, as they overtake the body. However, 
the drag remains small, and possibly tends to zero in the long time limit. The behavior of 
supercritical flows (Figure [T]) is quite different. A down-stream propagating solitary wave, with 
a small oscillatory tail, is generated and travels at constant velocity cq « —1/2. It remains 
a symmetric, depression wave at the location of the body, that does not generate any wave 
resistance. Unsurprisingly, this corresponds to the observations of the one-layer theory in the 
supercritical case. 

4.3.2 Critical Froude numbers 





(a) a = £2 = M = 0.1, <5 = 5/12, 7 = 0.99. 



(b) a = £2 = M = 0.1, S = 12/5, 7 = 0.99. 
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(c) Q = £2 = /i = 0.01, 5 = 5/12, 7 = 0.99. (d) a = £2 = m = 0.01, .5 = 12/5, 7 = 0.99. 

Figure 8: Dependence of the wave resistance coefficient on the Froude number, at time T = 10. 



As we have seen in the previous section, the solution of the fKdV equation, and therefore 



the wave resistance coefficient (jC.5|) . are small when the Froude number is away from its critical 
number Fr = 1. We present in Figure [8] the behavior of the wave resistance at time t = 10, 
depending on the Froude number, for different values of the depth ratio 6, and different values of 
the shallowness parameter /i. 

It appears that the maximum of the critical wave is obtained below the critical speed, at 
approximate value Fr « 1 — Cstt A4. = 1 — Cstt fj, ^^^g ■ The fact that the maximum peak is 
slightly subcritical has been obtained using the linear theory in [44| , and explained as a result of 
dispersion effects (without the nonlinear and dispersion terms, the peak is infinite and obtained 
at Fr = 1). This effect is therefore preserved in the nonlinear theory. 
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In Figure [HJ we compute the solution of (|4.16p . with zero initial data, in the critical case 
(Fr — 1). Again, we plot the flow, depending on space and time variables, as well as the wave 
resistance coefficient, and we set S G {5/12, 12/5}. It is known, since the work of Wu [5S], that the 
forced KdV equation can generate up-stream propagating solitary waves. Our simulations exhibit 
that behavior, and show that the depth ratio i5 now plays an essential role, as it determines the 
sign of the coefficient of nonlinearity. Indeed, even if up-stream propagating waves are generated 
for both of the values of S, these waves are of elevation if (5^ > 7, and of depression in the opposite 
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[ Wave resistance coefficient 










(a) 5 = 5/12. 



(b) 5 = 12/5. 



Figure 9: Critical flow: Fr = 1. a — €2 — — 0.1, 7 = 0.99. 



The continuous generation of up-stream propagating waves induce a periodic aspect to the 
wave resistance, that recalls the hysteretic behavior of the dead-water phenomenon recorded 
in [12] • However, we do not believe that these waves are the cause of the periodic behavior 
recorded in [32], based on the two following facts. 

i. The time period of the generation of waves AT is very large. A simple scaling argument 
shows that AT oc l//.t if c = and X,a,v cc /i, and the proportionality constant appears to 
be big (see also scaling arguments in [53 )• By Proposition 14.21 the KdV approximation is 
a valid model only up to times of order 0(1//!), so that nothing shows that the solutions 
of the full Euler equations (|2.4p follow this behavior. It is possible that the time-range of 
validity of the model is shorter than one time period of generation of up-stream propagating 
waves. 

ii. As a matter of fact, the solutions of the strongly nonlinear model p.3p do not exhibit such 
a phenomenon, as the generated internal wave stays in the trail of the body, and never 
overtakes the ship, even for long times. 

In Figure [H] page [THl we present the results of a simulation where the velocity of the ship is 
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affected by the generated wave resistance. A similar periodic behavior is exhibited, with a more 
reasonable time period (AT « 10). What is more, in this case, the up-stream propagating solitons 
vanish as soon as they overtake and slow down the body. This phenomenon corresponds to the 
observations of [12], but cannot be seen in Figure [HI This is why we believe that constant- velocity 
models do not accurately predict the hysteretic aspect of the dead-water phenomenon, and that 
only a constant-force model would be able to recover this behavior (see also discussion in 
section 1.3). 

5 Overview of results and discussion 

We have presented several nonlinear models describing the behavior of a system composed of two 
homogeneous fluids, when a body moves at the surface with constant velocity. Our asymptotic 
models are obtained through smallness assumptions on parameters of the system (Regime [1] and [H 
see page [TT|l . Contrarily to (linear) models existing in the literature, our models are valid for 
large (Regime [ij and moderate (Regime [2]) amplitude waves. They are justified by consistency 
results (Proposition l3.1|) or convergence results fPropositions 14.1] and I4.2( l . 

These models duly recover the key features of the dead-water phenomenon, described in details 
in [32], namely: 

i. transverse internal waves are generated at the rear of a body while moving at the surface; 

ii. the body suffers from a positive drag when an internal elevation wave is located at its stern; 

iii. this effect is strong only near critical Froude numbers; 

iv. the maximum peak of the drag is reached at slightly subcritical values. 

Moreover, we numerically computed our models using several ratios between the depth of the 
two layers of fluid. It appears that when the upper layer is thicker than the lower layer, the 
generated wave is able to reach larger amplitudes, and the wave resistance suffered by the body is 
stronger than in the opposite case (see figures [3] [Hand [5]). This aspect have never been expressed 
before, as far as we know. 

One important characteristic of the dead-water phenomenon, investigated in |42J, is the hys- 
teretic behavior of the system, and is not predicted by our models. This is due to the fact that our 
analysis assumes the velocity of the body to be constant, whereas a constant force was imposed 
to the body in their experimental setting. It would be of great interest to construct and justify 
models in the constant-force setting, in order to recover the complex behavior of the system in 
that configuration. 

Finally, our work is limited to the two-dimensional configuration, that is horizontal dimension 
d = 1. Extending the study to the three-dimensional configuration would allow to observe 
divergent waves generated by the moving body, as well as the transverse ones [441 159j . 

A Derivation of the Green-Naghdi type model 

In this section, we derive the strongly nonlinear model (j2.7|) . This model only requires the shallow 
water assumption 

and is a convenient intermediate step to construct the models of Regimes [T] and [2] used throughout 
the paper. The full Euler system is proved to be consistent with our model at order 0{fP), in 
Proposition lA.ll below. 
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Let us first plug the expansions of Lemma [2.61 into (|2.4p . and drop 0{ij?) terms. One can 
easily check that we obtain the approximate system 

{dt - Frdx)(2 + dx{h2dx^2) = fJ'dxT[h2 , 0]dxi^2 , 



(A.l) 



{dt - Frdx) {dx4'2 - idx^i) + (7 + S)dxC2 + jdx {\dxi^2? - 7|5.V'i|') 

fi(jidt - Ft: dx)dxU{dxi)i,dx4>2) + t2dx'R-2[dxi}i,dxil}2]^ 

dxC2 



^8^ 



Bo nv/l + Mei|5.C2P 



where we have used the following notations: 



T[h,b]V = -^dxih^dxV) + ^{dx(h\dxb)V)) - h^{dxb){dxV)) + h{dxbfV, 



M[Vi,V2_ 



{{dxh)V2 - dx{hV2)f ~ l{{dxh)Vr - dx{hV2)f 



■H{Vi,V2) = hi(dx{hiVi) + dx{h2V2)-hiidxVi-dx{hi+h2)Vi), 



r[hi,h2]Vi+r[h2My2~ l^dx{hldx{h2V2))-hidxh2dx{h2V2)), 



ni(Vi,V2) 

n2[Vi,V2] = iVidx(n{Vi,V2)) +M[VuV2]. 



Now, thanks to the fact that Ci is a forced parameter of our problem, the system reduces to two 
evolution equations for {(^2, v), with v the shear velocity defined by 



V = d. 



'x {(P2 



1 L = 8 



) = 0x^)2 -^H{%l;i,ip2)- 



From the last estimate of Lemma 12.61 one has immediately 

(A.2) V = dx^p2-ldxA-l^ldxnidxA,d.4'2)+Oi^i^). 

Now, one deduces from the first equation of (|A.ip . and (IA.2I) . the following relations 

hidxi^i + h2dx^J2 = aFrCi + 5,?/'2) 
(A.3) V = 0x^2 -idxA - t^ldxn{dxA,dxi^2) 

It follows that any linear operator defined above can be approximated as in the following example 

ft.2f — QfFr^i /iiu + 7aFrCi~ 



'H{dx1pl,dxTl'2) 



n 



hi + 7/12 hi + 7/12 

For the sake of readability, we do not precise the arguments in the following, and simply write 



n 



n 



TZi = Til 



h2V~aFT(^i hiv + ^aFT(^i 



hi + 7/12 



hi + 7/12 hi + 7/12 

Using (jA.3p . one can approximate dxi/Ji and dx'4'2 at order 0{^^) with 

{hi + -fh2)dx^i = ~h2V + aFrC + ^i (Ui - -/h2dxn) + 0(^1^) 
(A.4) {hi+jh2)dx^2 = hiv + jaFiC + n-f {TZi + hidxH) + Oifx^). 



hi + 7/12 
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Using these formulae, the system (lA.ll) becomes (dropping 0{ix^) terms) 

h2 



(5t-Fr9,)C2 + 9. 



hi + 7/12 



(A.5) <^ 



[dt - Frdoojv + (7 + S)dxC2 + 77^^= 



(/ii +7/12)2 
1 



-a 



dxC2 



where C and Q are defined by 

/l2 



£= 7 



/ii + 7/12 



{Ui + hidx-H) - r[/i2,o] 



Bo "Vx/i + Mei|5.C2P 

hiv + 7a Fr Ci 



(/^i + fe2)t^-(l~7)«FrCi ^ ^ 



(/ii +7/^2)2 



/Ji + 7^2 
ft,2W~aFrCi hiv + "faFT(^i 

^1+7^2 ' /ll+7/l2 



((fef - 7fei)f + 7(/ii + h2)a Ft (i ^ 



(/^i+7/i2)2 

It is now simple to check that the system (jA.5|) is exactly the system p.7p . 

This model is justified by the following proposition: 

Proposition A.l. The full Euler system (j2.4l) is consistent with (jA.Sp . precision ©(/i^) on 
[0,r], TOi/i T > 0. 

Proo/. Let [/ = ((1,(2, "01, "02) be a strong solution of bounded in W'^'°°i[0, T];H''+*'>) with 

s > 1 and to > 9/2 + 5, and such that (j2.5p is satisfied with (i{t, x) = (i{x). Using Lemma [2Jl 
it is clear that (Ci, C2, 9^-01, 5^-02) satisfies (|XT|) . up to Ri = (ri, ra, rg, r4)^ e L°°([0, T]; 
satisfying (for i = 1 . . . 4) 



l^lwi'~H=+*o 



Now, we set w = dxi^2 - 7-ff(?/'i, ^-2), and v satisfies (|X2|) . up to R2 G iyi'°°([0, T]; iJ'^+S), with 
(again, thanks to Lemma 



3+tO 



Therefore, since TZi and S^jH involve two spatial derivatives, one has 

{hi + jh2)dxipi = -h2V + aFrC + fJ-Rs, 
{hi + -fh2)dxip2 — hiv + jaFr ( + ^R4, 

with i?3, Ri e L°°([0,T];i7*+3). Consequently, using the fact that H%R) is an algebra for 
s > 1/2 and that (^3)) is satisfied, we deduce that (|X4| is satisfied, up to R5 e L°°{[0, T];H^+^), 
with 

1 



U?5 



i?4 



Finally, when plugging {v,C,2) into (lA.ip . the residuals are clearly bounded by /i^Co, uniformly 
in L°°{[0, T]; iJ*), and the proposition is proved. □ 
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B Proof of Proposition 14.21 

The sketch of the proof is the following. We will first prove the convergence between 



defined by (|4.8p - (|4.11l) . and the solutions of the symmetric Boussinesq-type system (|4.3p . This is 
obtained thanks to a consistency result, together with energy estimates. The convergence towards 
solutions of the full Euler system (|2.4p follows then immediately from Proposition 14.11 Finally, 
the proposition is completed when remarking that the corrector term Ui obeys to a sublinear 
growth. This result is a consequence of the following Lemma, that proceeds from Propositions 3.2 
and 3.5 of 



Lemma B.l. Let u he the solution of 
{dt + cdx)u = g{vi,V2), 



(B.l) 



with Vi e {1,2}, 



with ci 7^ C2, Vi, V2 G -ff^(K), s > 1/2, and g is a bilinear mapping defined on R and with values 
in R. Then one has the following estimates: 

I. Ifc^ci , then ^hin =0. 
a. Ifc^ci^ C2, then ■^\u{t, ^[/^.(r) = ^'(1). 

Moreover, if there exists a > 1/2 such that Vi{l + x^Y" , and ^2(1 + G H^{R), then one has 
the better estimate 



< Coy,{i + x')X.M{i + -')XHm^ 



with Cq = C(c, ci , C2 ) . 

We can now proceed with the proof. 

Step 1: Well-posedness of Uapp{t, x) . The global well-posedness of the forced Korteweg-de Vries 
equation is given by Bona and Zhang in [5]. 

The proof relies on an a priori estimate, that we recall here, as it will be useful for the following 
arguments. Let u be a solution of 

dtu + cdxU + XudxU + i^d^u = b{x). 

As we multiply the equation by A^'' u (with s' > 3/2), and integrate with respect to the space 
variable, one obtains 



A / A** {udxu)A'' u dx 



(A" 6) (A" u) dx 



Thanks to the Kato-Ponce Lemma, we estimate the right-hand side as follows: 



A^ {udxu)A'^ u dx 



< 



[A' ,u]dxu{^' u) dx 



^1 I ^^l^l 



and the Cauchy-Schwarz inequality leads to 



1 d 



2dV 1^ 



\hA^\h^ 
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It follows 



sup \U 

te[o.T] 



Hs' <exp^A J \d^u{s,-)\^^ dsj |^,, +^Hi/= 



In the particular case of (|4.12D . since a = 0{^), 62 = and ri± — 5(C2 ± ^^''^) \t=o ■ 

it is straightforward that there exists T = T(7 + (5, ^^qrj) > such that for all s' > 3/2, 

(B.2) sup^ ^\u\^^, <Co(|^Ci|^.',|^L.o \h^',T). 



te[o,T/H 



We can then exhibit Ui. Indeed, (j4.11l) can be written in a simplified form as 
(B.3) {dt + Cidx)e^ ■ SoUi = ^ fjjk (r, t , a;) + d^gij (t, t , , 

with 

fijk{T,t,x) EE a^jkUj{T,x - Cjt)d.j:UkiT,x - Ckt) and gij{T,t,x) = l3ijdluj{T, x ~ Cjt). 
Thanks to estimate (|B.2|) . one has for s' > 1/2, 

|/»jfeLoo([o,T]x[0,T/e);ff»' + i) + \9v\l'=°([o,T];H''') - 1 H='+2 ' 1^ \t=o |fl'='+2,T')- 

with Co = Co(::^,7 + 5). Hence, for s' > 1/2, one can set 

e± ■ SoUi{T,t,x) ^ ^ / /yfc(T, s,x + Ci(s-t)) ds 

+ X! ^ (Sij ("T: 2; - Cjt) - (r, X - Cii)) 



(B.4) = ^ W^''+^V'^. 

One checks immediately that f7i e L°°([0,r] x [0, T/^); i?"') and iJi = 0. 

Ste'p 2: Estimate on Ui . Let us estimate each term of the decomposition (|B.4p . Thanks to the 
uniform estimates of Step 1. above, one has gij S L°°{[0,T]; H'^ ), so that it follows immediately 

Vj^i, |^''L~(R+x[0,T];H=') ^ '^o(|^Cl|ff.'+2,|V'L=o |/fs' + 2,7'). 

Moreover, for j 7^ i, we remark that fijj can be written as 

SO that has actually the same form as F*-', and can be treated in the same way. Since 
f,,j e L-([0,r] X [0,T/e);i/^'+i), C/^^'^" is bounded in H^' by Co(|£Ci|h.', 1^^ L=o 1,/.'.?^): for 
ah 

Finally, for all {j,k) 7^ with j 7^ fc, t/*^'^ satisfies the hypothesis of Lemma [B.ll with 
fijk — g{uj,dxUk). Therefore, we deduce 

(B.5) l^iLoo(R+x[o,T];if=') ^ ^C'o(|^Ci|^.'+2, |^^L=o \h^'+2,T). 
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As for the second estimate of the proposition, let us first remark that the estimates of V^^ 
and W"^^ are time-independent, and in agreement with the improved estimate. Therefore, the 
only remaining terms we have to control are C/*^'^ with j ^ k. Of course, we will use the second 
case of Lemma FB.li but we have to check first that for every r e M+, the initial data Wj(t, 0, x) 
and dxUk^T, 0, x) are localized in space, that is 

VreM+, 1(1 +a::^)u±(T,0,a;)|^,, + \{1 + x'^)dxU±{T,0,x)\^^, < oo. 

This property is true at t = (by hypothesis of the proposition), and is propagated to r > 0, using 
the fact that u± (r, x± ) satisfies the KdV equation (|4.10p . This propagation of the localization 
in space has been proved for by Schneider and Wayne in [SH Lemma 6.4]. We do not recall the 
proof here, and make direct use of the statement: 

Lemma B.2. If {1 + x^)V € i/* then there exists Ci,Ci > such that 
|(l + a;2)u±(r,0,x)|^,,+, < Ci |(1 -f 2;2)u± |_^^ < Ci|(l + x^)^ 

This Lemma, together with the second estimate of Lemma IB.li allows to control uniformly 
jjijk g ^oo^]^+ X [o,r];iJ«'). Every term of the decomposition (iRi)) has been controlled, and 
one has the following estimate: 

(B.6) |t^iL~(R+x[o.T];H-') - '^o(|^Ci|^.'+2, |^L=o Ih='+2, |(1 + a;^)^L=o In^'+^^^y 

Step 3: Consistency results. We first prove the the consistency of the Boussinesq-type system (|4.3p 
with our approximation. The precision is 0{y?/'^) in the general case, and 0[y?) in the second 
case. Here and in the following, we set a = €2 = in order to simplify the notations. The general 
case a = 0{fi), 62 = is obtained with slight obvious modifications. 

Plugging Uapp{t,x) into (|4.6p . we see from (|4.7p - (|4.1ip that the only remaining term we have 
to control is ^^R{^t,t,x), with 

R = drUi + J:iiUo)dxUi + ^iiUi)dxUo + SiiUo)dtUi + SiiUi)dtUo 
-^2dlUi - S^dldtUi -f /iSi(C/i)a,C/i + tiSi{Ui)dxUi, 

where Uo{fj,t,t,x) = u+{^t^x — c+t)ei + u^{fit,x — c_i)e_. We bound each term of the 
right-hand side in the Sobolev iJ^-norm, with s > 1/2 as in the proposition. 
The estimate (jB.Sp . with s' = s + 3, leads to 

\^2^lU^\^^ < Vt Coi\-^Ci\H.+.^Vl^o 
Then, from (|4.1ip . one has 

• dtSoUi = -CiGi ■ dxSoUi + fi 
with /, e L-([0,r] X [0,T//i);i7^+2), and |/,|^,^, < Co(|£Ci|j^.+5, I^^L^o \h^+..T), so that 

One obtains in the same way the desired estimates for Si(C/o)9a;t/i, T,i{Ui)dxUo, Si{Uo)dtUi, 
SiiUi)dtUo, Ei(C/i)a,[/i and Si{Ui)dxUi. 
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Finally, in order to estimate drlli, we differentiate (|4.1ip with respect to r. Since Ui satis- 
fies (|4.10p . one has drUi G L°°{[0, T]; We are on the framework of Lemma fB.li so that we 
can obtain, as for that drUi G L°°([0, T]; with 

Hence, the residual fi^R is uniformly bounded in L°°([0, r//i); iJ"), and 

which gives the consistency at order ©(/i"^/^). 

The consistency at order 0{p'^), for initial data sufficiently decreasing in space, is obtained in 
the same way, using the second estimate (jB.6|) . 

Step 4- Convergence results. The convergence is deduced from the consistency, as in the proof of 
Proposition HTT] Indeed, when setting Rf^ = Uapp — V = Uo{fJ.t, t, x) + iJ.Ui{fj.t, t, x) ~ V{t, x), one 
can check that i?'^ satisfies the following equation: 

(B.7) {So - fiS^dl + tiSiiU,pp)dtR'' + (So - ^^25^ + ^Ei(C/app))a,i?'' = ^'^V + M + M^, 

with A = 5f5i([/app)i?^ - Si{Ri')dtV, B = d^^i{Uapp)Rt' - i:i{Ri')d^V and the function / 
proved to be uniformly bounded in by the consistency result. We define the energy as 

= i(5oA^■i^^ A^i?'^) + |(52a,A^i?^,a,A^i?^) + |(5i(C/app)A'^i?^ A^i?^), 

as in the proof of Proposition 14.11 The same calculations lead to the following differential in- 
equality: 

^E.iR'^) < Co^iE,iRn + CoM'/2(i^.(i?^))^ 
at 

with Co = C'o(|g^Ci|//s+5j |^lt=o l^f^+si^^)- Then, Gronwall-Bihari's theorem allows to obtain 
{EsiR^'))^/^ < CoAi^/^(e^o''* - 1), and finally for ^it < T(:^,7 5), 

|C/app-l^|^. <Co(i?.(i?^))i/2 <M'/'t Co(|^Ci|^.+.,|^L.„ 
The first estimate of the proposition is now a direct consequence of (jB.sp . and 

|^Ci(a;)|jj,+5 + |^^(^,a;)|^.+5 < \U{t,x)\^^^,^, 

thanks to an appropriate estimate of the operator i/(-0i,-02); see [23l Proposition 2.7]. 

The second estimate of the proposition follows in the same way, using (jB.6p and the consistency 
at order 0{iJ?). 



C Wave resistance and the dead-water phenomenon 

Following Lamb |32| and Kostyukov 31 , we assume that the drag experienced by ships is mostly 
due to the wave (making) resistance. This section is devoted to the analysis of this wave resis- 
tance. We first work with the variables with dimension^ and deduce an explicit formula for the 
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wave resistance Rw, depending on the flow, given as a solution of the full Euler equation (|2.ip . 
Accordingly with the nondimensionalization performed in Section 12.21 we introduce the dimen- 
sionless version of the wave resistance, that we call wave resistance coefficient and denote Cw- 
Finally, we derive simple approximations, for the two regimes considered throughout the paper. 

The wave resistance acknowledges the energy required from the body to push the water away, 
and is defined by (see [151 1331 [ST] and references therein) 

Rw ^ [ P(-e, •n)d5, 

where Fship is the exterior domain of the ship, P is the pressure, is the horizontal unit vector 
and n the normal unit vector exterior to the ship. The pressure being constant (= Poo) on the 
non-submerged part of the ship, one has 

Rw^ f (P-Poo)(-e, -n) d5 = / -Poo)(-9.Ci) dx 

= - / ^'Li+ci dx. 
Jr 

As a solution of the Bernoulli equation in (j2.ip . the pressure P is given in the upper domain by 
P(x, z) 1 

— = -dt(j)i{x,z) - -|Va;,z(^i(a;, z)p - gz, 

Pi 2 

so that the wave resistance satisfies 

Rw ^ Pi J g{di + Ci)9xCi - {dxdtcjii + ]^d,^ {d,j:4>if + ]^dx {dz4>if^ Udi+ci Ci da;. 

= / PH'[C2,</'l](dl+Cl)Cl dx. 
jR 



Let us now construct the dimensionless version of this formula. Using the same change of 
variables as in Section [2.21 it is straightforward to obtain the dimensionless wave resistance (that 
we call wave resistance coefficient, and denote Cw)'- 

(C.l) Cw = Rw ^ I Fw[k]{i)C,i{i) df, with 

-Fwi^i] = {di~Fidi)di4>i + jdi(^(^di^iy + ^(^di4>iy 

Again, we omit the tildes in the following, for readability. Let us remark that by definition of the 
Dirichlet-Neumann operator Gi, and using (j2.4l) . one has 

(C.2) - Fw[<l^i] = [dt - Fra,)9,0i + ({d^^if + (^-aFra^Ci + ei^Ci^.^i] J • 



For practical purposes, we use approximations to compute the wave resistance coefficient, 
corresponding to the leading order term in the asymptotic expansion of Cw, for the different 
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regimes at stake. These formulae use the variables of Sections [3] and HI i.e. the surface and 
interface deviations and C2 and the shear velocity, defined by 

V = 9a; (02 -70lL=.3cJ- 

First of all, we use the following estimate, justified in [53]: 

cj)i[x,z) = i'lix) + 0{y). 



When combined with the first equation of (jA.4[) : 

— h'2V ~1~ Fr 1 

dx^i = r — ; — r +0{^JL), with /ii = l + eiCi, and /i2 = T+e2C2, 

hi + 7/12 

one has immediately that (jC.2p simplifies into 

-h2V + a Fr 



hi + 7/1: 



2 



_,e2. /^ (-fe2i^ + aFrCi)V(a£Ci)M(fei+ 7/^2) Fr +62(^2^ -«FrCi))' ^ , , 

Now, 9t/i2, and can be written using only C,i, ^2, and their spatial derivatives, since they 
satisfy system (lA.Sp up to C(/i^). Therefore, the leading order term of the wave resistance 
coefficient Cw can be deduced from the knowledge of the solution {C,2,v). We do not write 
explicitly this expression, as the models used in our simulations benefit from extra smallness 
assumptions, and simpler formulae are deduced in these cases. 



Case of Regime [T] ; 



M « 1 ; a = ^ = o{^l) , 1-7 = 0{^i). 
(2 



The first immediate simplification in this regime is 

/ll+7/l2 = /ii + /l2 + 0{^i) = \ + \ + 0{fi). 



Then, we use that {C2,v) satisfies system p.ip up to 0{fi^) to deduceH (with hi = I + j — /12) 

dt{h2v) = FYdx{h2v) + (1 + d)h2dxC2 + £2 (^vd.. {hih2v) + yS, {hi - /la)^')^ + 0{^l) 

= Yrdx{h2v) + (1 + 5)h2d,h2 + e2dx{h2v) + ^-L^d^i^k^), 



Finally, using these approximations into (IC.ip and (|C.3|) . one has 

Cw = / f (1 + 5)/^29xC2 + e2dx{h2v) + e2-^c'.(/i2^'')) (1(2^) + ©(/x^) 

^as we use Bo^^ = fj? in our simulations, we do not take here into account the surface tension term. 
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This is the formula used in Figures [5HS1 in Section 13.21 
Case of Regime [2] : 

< 1 ; £2 = , a ^ 0{ji). 

Most of the terms of (jC.3l) are now of size ©(/i). The first order of system (j3.3l) leads immediately 
5 to 

dt{h2v)~Yid^{h2v) = -2±ld.X2+0{y), 



so that we obtain simply 

(C.5) Cw = - U2^ClAx + o{^i). 

Jr da; 

This formula is used in Figures ISHHl in Section 14.31 

10 Remark C.l. As we can see in (|C.4p and (jC.Sp . the wave resistance coefficient will he small in 
the two following cases 

i. The internal wave is small, or is not located below the ship, 

ii. The internal wave is symmetric and centered at the location of the ship. 

The first case is obvious, and the last case reflects the fact that the integrands of (jC.4|) and (IC.SP 

15 are odd if C2, ^ cind are even. 

On the contrary, the ship will encounter a strong positive wave resistance if the functions (2 
and h2V are decreasing at the location of the body. This is the case when an internal wave of 
elevation is located just behind the ship. The dead-water effect is then explained in that way: the 
ship, traveling in a density- stratified water, generates internal waves of elevation in its trail, and 

20 therefore suffers from a severe wave resistance. 



D The Numerical schemes 

We present in this section the numerical method used to obtains the figures of this study. For 
all of the simulations, we use a scheme based one the Crank-Nicholson method, and take care of 
the nonlinearities using a predictive step. This method has been introduced in [U [3] , and used 
in the water wave framework in [15l [24l [22] . We give the general directions of such a method in 
the following, and then present the exact schemes we used for each of the models. 

Time discretization. Denoting At the time step of the scheme, we approximate a function u{x, t) 
at time t = uAt by u{x, uAt) = Then, we approximate the time derivative by 

dtu : 



At 

and any linear function of u by 



^^^^ _ i^(u"+i) + F(m") _ 



2 V 2 

Now we deal with the nonlinearities using a predictive step, u"^^ , defined by 



u ' 2 + u 2 
= w • 
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The first half-step, W2 , is computed through a simple explicit scheme. 

With the help of the predictive step, the discretization of a quadratic nonlinearity F(u, v) 
(F being linear with respect to both variables) is then a linear combination of the two possible 
discretizations, namely 

;"+^ j +(l-a)F(^ ,i;"+2j . 

The parameter a can be chosen so that natural quantities are conserved. 

Space discretization. With Ax the space step of the scheme, the functions are discretized spatially 
with a central difference. It is also useful to introduce a Lax-type averaging, by 

u{x,t) « M{f3)u with (M(/3)u"), = (1 - /3)< + ^(u^, + u^_,). 

The spatial derivatives are given by 

2Ai = (^i")», d.,u^ — = {D2u)„ ... 

We use periodic boundary conditions. 

As we see in the following section, the parameters a and /3 can be set so that natural quantities 
are conserved. In the case of a Korteweg-de Vries equation, when choosing carefully the param- 
eters, one can obtain a numerical scheme that preserves the discrete ^^-norm, in agreement with 
the solutions of the KdV equation having constant L^-norms. As for the more complicated case 
of system p.3p . we set the parameters by analogy with the simple Korteweg-de- Vries equation. 



F{u,v) w aF 



D.l The forced Korteweg-de- Vries equation 

We want to solve numerically the following generic forced KdV equation 

(D.l) dtu + cdxU + XudxU + lyd^u = f {x) . 

When there is no forcing term (/ = 0), it is known that the KdV equation preserves the energy 

I 1 2 

E = |m|^2- The following scheme has been presented and studied (without the forcing term) 
in |24j . and has the property to preserve the discrete energy when / = 0. 




Proposition D.l. If f — 0, then the scheme (|D.2p preserves the discrete P-norm: 

i i 

The proof is given in [23] Theorem 2] . 
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D.2 The fully nonlinear model of Regime [T] 



The system we deal with now is (I3.3p . that we can write under the compact form 
(D.3) 

idt-FTd,)iU-fini[U]) + iiA[U] + eiBix))U) + Me2a,(7^2[C/] • C/) = b^x) + ^dl{T[U]) . 

DO 



withC/ = {C2, w)^,b{x) = -aFr/(l + ,5) (Ci(x),0)^, and 



+ S e2jf^f^^wj h \ U Qi{x)J 



It is convenient to denote, with hi — 1 + ei^i — e2C2 and h2 — j + €2(2, 



= T — ^ — 5K2] = £277 — \ — rri, and t[C2] = 



10 Advised by the above work on the KdV equation, we use the following time discretization: 



/[C]w«/[r+^] 7, : 



1 + w 



--dx 



Qn+l 



dt{Si[h]w) ^ Si[h"^ 



At 



/■n+1 _ /-n 



with 5i[/i,u;]C = jdl{C{l + l/6-2h)w)-2e2{dxC){dxh)w, 



W S2[H]w « -w"+2 52[H"+2] ( 



2 + ^"+1 
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Finally, after the space discretization, this leads to the following scheme (with ft," = | + e2C" 
(D.4) 

.„+i _ , +1 j^(-n\ Fr J / / / 1 C"+^ + C" \ 



+ fM(l)/[/i"+3] ^ + /[ft"+3]Af(l/2)(u;" + u;"+i; 

(D.5) 



+ (7 + -5) 



'1- 



Fr(5 / /7/;"+i + 
^L*! M(l)g[/i"+^]7«"+3 +.g[/i"+5]u;"+H/(l/2)(u;" + u;"+i) 

' /-n+l I Ar; 




n+i|2 



1 + /iell^xC 



D.3 Validation of the method 

In order to validate the proposed schemes, we use known explicit solutions of the forced Korteweg- 
de Vries equation (|4.13p . The first test function we use is the travelling wave 

(D.6) uq = — sech^(fc(a; — Cswt)), 



with Csw = c + Avk^, k = ^y—X/{12l'), corresponding to the classical case where there is no 
forcing. 

The second test function is the steady solution 
(D.7) Ml = ±sech^(fca;), 

with fc^ — and the corresponding forcing 



f{x) — (c + — ^ sech^(fcx). 



These two functions are exact solutions of the forced KdV equations 

dtu + cdxU + XudxU + vdtu — —f(x). 

ax 

In the absence of non-trivial solutions of system p.3p . we use the same functions as reference. 
After an appropriate change of variables, uq and ui will satisfy system (13.31) up to small terms. 
Therefore, adding the corresponding forcing term to the equations, uq and ui are exact solutions 
of the modified system, and we are able to compare the results of our numerical schemes with 
the theoretical solution. 

These comparisons are presented in Tables [5] and |31 The error are givens in term of the 
normalized P norm. The results exhibit a convergence behavior of order ©((Ax)^ -f (Ai)^). 
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Ax — At L T KdV scheme fully nonlinear scheme 

~~0A 20 To 8.1530 lO^^ 4.9498 lO-"* 

0.05 20 10 2.0393 lO""^ 1.5154 lO^'' 

0.01 20 10 8.1604 10-6 1.8363 10"^ 



Table 2: Numerical errors of the KdV and Green-Naghdi schemes for the initial data (jD.6|) . We 
choose the parameters fi = €2 = a = 0.1, 7 = 0.9, 6 = 5/12, Fr — 1.1, Bo — 100. 



Aa; = At L T KdV scheme fully nonlinear scheme 

~~0A 20 To 4.7397 lO"* 1.7844 lO""* 

0.05 20 10 1.1850 10--* 4.4557 10"^ 

0.01 20 10 4.7409 lO^^ 8.1604 10"^ 



Table 3: Numerical errors of the KdV and Green-Naghdi schemes for the initial data (|D.7p . We 
choose the parameters — €2 — a — 0.1, 7 — 0.9, 5 — 5/12, Fr — 1.1, Bo — 100. 

The test functions ([HB and (IDJ)) are very smooth, and our problem can lead to rapid 
variation in the solution, as we can see especially in figures [5] and [71 In order to check the 
accuracy of our method in such cases, we compute the KdV scheme with different space and 
time steps, and compare the outcome with the result of a more refined grid. More precisely, we 
compute the scheme for Aa; = At = 0.1, 0.05, 0.01, and measure the difference with a reference 
solution obtained with Ax — At — 0.001. The error is then discrete P norm of the difference, 
normalized with the P norm of the reference solution. 

These comparisons, using the settings of figures [5] and [3 are presented in Table S) The 
normalized error can be relatively large; this is due to the fact that the reference solution is very 
small (of order 10^^). However, one clearly sees that the schemes are stable, and converge fast 
when space and time steps become small. Moreover, the precision for Ax = At = 0.01 is sufficient 
to validate our results. 



Ax = At 


Figure [5] (a) 


Figure El (b) 


Figure [7] (a) 


Figure [71(b) 


0.1 


7.2320 10-1 


1.2158 


2.4538 10-1 


2.0070 10-1 


0.05 


1.9339 10-1 


5.1843 10-1 


9.6817 10-2 


1.0586 10-1 


0.01 


9.3586 10-3 


5.2038 10-2 


1.9121 10-2 


6.0949 10-2 



Table 4: Convergence of the KdV scheme, with the settings of figures [6] and [T] 
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